source: src/Actions/FragmentationAction/FragmentationAutomationAction.cpp

Candidate_v1.6.1
Last change on this file was a6c11a, checked in by Frederik Heber <frederik.heber@…>, 7 years ago

All failures in FragmentationAutomation are explained by log statements.

  • Property mode set to 100644
File size: 21.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 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * FragmentationAutomationAction.cpp
26 *
27 * Created on: May 18, 2012
28 * Author: heber
29 */
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36#include <boost/archive/text_iarchive.hpp>
37// boost asio needs specific operator new
38#include <boost/asio.hpp>
39
40//#include "CodePatterns/MemDebug.hpp"
41
42//// include headers that implement a archive in simple text format
43#include <boost/archive/text_oarchive.hpp>
44#include <boost/archive/text_iarchive.hpp>
45
46//
47//#include <boost/mpl/remove.hpp>
48//#include <boost/lambda/lambda.hpp>
49
50//#include <iostream>
51
52#include "CodePatterns/Assert.hpp"
53#include "CodePatterns/Info.hpp"
54#include "CodePatterns/Log.hpp"
55
56#ifdef HAVE_JOBMARKET
57#include "JobMarket/Jobs/FragmentJob.hpp"
58#else
59#include "Jobs/JobMarket/FragmentJob.hpp"
60#endif
61
62#include "Fragmentation/Automation/FragmentJobQueue.hpp"
63#ifdef HAVE_JOBMARKET
64#include "Fragmentation/Automation/MPQCFragmentController.hpp"
65#else
66#include "Fragmentation/Automation/MPQCCommandFragmentController.hpp"
67#endif
68#include "Fragmentation/Exporters/ExportGraph_ToJobs.hpp"
69#include "Fragmentation/Summation/Containers/FragmentationChargeDensity.hpp"
70#include "Fragmentation/Summation/Containers/FragmentationLongRangeResults.hpp"
71#include "Fragmentation/Summation/Containers/FragmentationResultContainer.hpp"
72#include "Fragmentation/Summation/Containers/FragmentationShortRangeResults.hpp"
73#include "Fragmentation/Summation/Containers/MPQCData.hpp"
74#include "Fragmentation/KeySetsContainer.hpp"
75#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
76#include "Fragmentation/Automation/VMGDebugGridFragmentController.hpp"
77#include "Fragmentation/Automation/VMGFragmentController.hpp"
78#include "Fragmentation/Summation/Containers/GridDownsampler.hpp"
79#include "Fragmentation/Summation/Containers/VMGData.hpp"
80#include "Fragmentation/Summation/Containers/VMGDataFused.hpp"
81#include "Fragmentation/Summation/Containers/VMGDataMap.hpp"
82#include "Fragmentation/Summation/Containers/VMGData_printKeyNames.hpp"
83#endif
84#include "World.hpp"
85
86#include <boost/bind.hpp>
87#include <boost/thread.hpp>
88#include <fstream>
89#include <iostream>
90#include <string>
91#include <vector>
92
93#include <boost/mpl/for_each.hpp>
94
95#include "Actions/FragmentationAction/FragmentationAutomationAction.hpp"
96
97using namespace MoleCuilder;
98
99// and construct the stuff
100#include "FragmentationAutomationAction.def"
101#include "Action_impl_pre.hpp"
102/** =========== define the function ====================== */
103
104class controller_AddOn;
105
106// needs to be defined for using the FragmentController
107controller_AddOn *getAddOn()
108{
109 return NULL;
110}
111
112static void updateSteps(Process &p, const size_t step, const size_t total)
113{
114 LOG(1, "There are " << step << " steps out of " << total << " done.");
115 p.setCurrStep(step);
116}
117
118template <class ResultClass>
119std::pair<double, double> getResiduals(const ResultClass &_result)
120{ return std::pair<double, double>(0.,0.); }
121
122std::pair<double, double> getResiduals(const MPQCData &_result)
123{ return std::make_pair(_result.accuracy, _result.desired_accuracy); }
124
125std::pair<double, double> getResiduals(const VMGData &_result)
126{
127 // sensibility check on absolute residual
128 if (_result.residual > 1e-5)
129 ELOG(2, "Encountered absolute residual greater than 1e-5: " << _result.residual);
130 // take smaller value of the two as just one of them needs to be less than precision
131 const double residual = std::min(_result.relative_residual, _result.residual);
132 return std::make_pair(residual, _result.precision);
133}
134
135template <class ResultClass>
136bool checkResults(const typename std::map<JobId_t, ResultClass> &_results)
137{
138 bool result = true;
139 for (typename std::map<JobId_t, ResultClass>::const_iterator iter = _results.begin();
140 iter != _results.end(); ++iter) {
141 std::pair<double, double> residuals = getResiduals(iter->second);
142 if (residuals.second != 0.) {
143 if (residuals.first >= residuals.second) {
144 ELOG(1, "Fragment Job " << iter->first << " converged to " << residuals.first
145 << " instead of " << residuals.second);
146 result = false;
147 }
148 } else {
149 LOG(4, "DEBUG: Not checking accuracy because desired precision not given.");
150 }
151 }
152 return result;
153}
154
155
156#ifdef HAVE_JOBMARKET
157/** This performs a direct sum calculation for the nuclei-nuclei interaction energy
158 * and forces for each fragment contained in \a shortrangedata.
159 *
160 * Results are (re)placed in \a _longrangedata.
161 *
162 * \warning This is only valid for fully open boundary conditions.
163 *
164 * \param _shortrangedata containing information on position and charges per fragment
165 * \param _longrangedata containing energy and forces on return
166 */
167/* Is not used currently, commented out to disable warning message.
168static void calculateNucleiNucleiLongRangeContribution(
169 const std::map<JobId_t, MPQCData> &_shortrangedata,
170 std::map<JobId_t, VMGData>& _longrangedata
171 )
172{
173 ASSERT( _shortrangedata.size() == _longrangedata.size(),
174 "calculateNucleiNucleiLongRangeContribution() - shortrange and longrange data have inequal size.");
175 std::map<JobId_t, VMGData>::iterator longrangeiter = _longrangedata.begin();
176 std::map<JobId_t, MPQCData>::const_iterator iter = _shortrangedata.begin();
177 for (;iter != _shortrangedata.end(); ++iter, ++longrangeiter) {
178 const MPQCData &data = iter->second;
179 VMGData &longrange_data = longrangeiter->second;
180 // set long-range contributions to zero
181 longrange_data.nuclei_long = 0.;
182 longrange_data.forces.clear();
183 longrange_data.forces.resize(data.positions.size(), FragmentForces::force_t(3,0.));
184 longrange_data.hasForces = true;
185 // go through positions and evaluate sum naively
186 ASSERT( data.positions.size() == data.charges.size(),
187 "calculateNucleiNucleiLongRangeContribution() - positions and charges differ in size.");
188 std::vector< std::vector<double> >::const_iterator positer = data.positions.begin();
189 std::vector< double >::const_iterator chargeiter = data.charges.begin();
190 FragmentForces::iterator forceiter = longrange_data.forces.begin();
191 for (; (positer != data.positions.end()) && (chargeiter != data.charges.end());
192 ++positer, ++chargeiter, ++forceiter) {
193 ASSERT( positer->size() == NDIM,
194 "calculateNucleiNucleiLongRangeContribution() - position "
195 +toString(std::distance(data.positions.begin(), positer))+" has not 3 components.");
196 const Vector position((*positer)[0], (*positer)[1], (*positer)[2]);
197
198 std::vector< std::vector<double> >::const_iterator otherpositer = data.positions.begin();
199 std::vector< double >::const_iterator otherchargeiter = data.charges.begin();
200 for (; (otherpositer != data.positions.end()) && (otherchargeiter != data.charges.end());
201 ++otherpositer, ++otherchargeiter) {
202 if ((positer == otherpositer) || (chargeiter == otherchargeiter))
203 continue;
204 ASSERT( otherpositer->size() == NDIM,
205 "calculateNucleiNucleiLongRangeContribution() - other position "
206 +toString(std::distance(data.positions.begin(), otherpositer))+" has not 3 components.");
207 const Vector otherposition((*otherpositer)[0], (*otherpositer)[1], (*otherpositer)[2]);
208 const Vector distance = position - otherposition;
209 const double invsqrdist = 1./distance.Norm();
210 const double factor = (*otherchargeiter)*invsqrdist*invsqrdist*invsqrdist;
211 (*forceiter)[0] = factor*distance[0];
212 (*forceiter)[1] = factor*distance[1];
213 (*forceiter)[2] = factor*distance[2];
214 longrange_data.nuclei_long += 0.5*(*chargeiter)*(*otherchargeiter)*invsqrdist;
215 }
216 }
217 }
218}
219*/
220#endif
221
222ActionState::ptr FragmentationFragmentationAutomationAction::performCall() {
223 boost::asio::io_service io_service;
224
225 // TODO: Have io_service run in second thread and merge with current again eventually
226
227 FragmentationResultContainer &container =
228 FragmentationResultContainer::getInstance();
229 const KeySetsContainer& keysets = FragmentJobQueue::getInstance().getKeySets();
230 const KeySetsContainer& forcekeysets = FragmentJobQueue::getInstance().getFullKeySets();
231 const FragmentJobQueue::edges_per_fragment_t &edges_per_fragment =
232 FragmentJobQueue::getInstance().getEdgesPerFragment();
233
234 size_t Exitflag = 0;
235 std::map<JobId_t, MPQCData> shortrangedata;
236#ifdef HAVE_JOBMARKET
237 {
238 const size_t NumberJobs = FragmentJobQueue::getInstance().size();
239 MPQCFragmentController mpqccontroller(io_service);
240 mpqccontroller.setHost(params.host.get());
241 mpqccontroller.setPort(params.port.get());
242 // Phase One: obtain ids
243 mpqccontroller.requestIds(NumberJobs);
244 if (mpqccontroller.getExitflag() != 0) {
245 ELOG(1, "Could not request all ids from Server.");
246 return Action::failure;
247 }
248
249 // Phase Two: add MPQCJobs and send
250 const bool AddJobsStatus = mpqccontroller.addJobsFromQueue(
251 params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
252 params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly
253 );
254 if (AddJobsStatus)
255 LOG(1, "INFO: Added jobs from FragmentJobsQueue.");
256 else {
257 ELOG(1, "Adding jobs failed.");
258 return Action::failure;
259 }
260 mpqccontroller.run();
261
262 // Phase Three: calculate result
263 setMaxSteps(NumberJobs);
264 mpqccontroller.setUpdateHandler(
265 boost::bind(&updateSteps, boost::ref(*this), _1, _2)
266 );
267 start();
268 boost::thread wait_thread(
269 boost::bind(&MPQCFragmentController::waitforResults, boost::ref(mpqccontroller), boost::cref(NumberJobs))
270 );
271 wait_thread.join();
272 stop();
273 if (mpqccontroller.getExitflag() != 0) {
274 ELOG(1, "Could not obtain all results from Server.");
275 return Action::failure;
276 }
277
278 mpqccontroller.getResults(shortrangedata);
279 if (!checkResults(shortrangedata)) {
280 ELOG(1, "At least one of the fragments' energy could not be properly minimized.");
281 return Action::failure;
282 }
283 if (mpqccontroller.getExitflag() != 0) {
284 ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag.");
285 return Action::failure;
286 }
287
288 Exitflag += mpqccontroller.getExitflag();
289 }
290#else
291 {
292 const size_t NumberJobs = FragmentJobQueue::getInstance().size();
293 MPQCCommandFragmentController mpqccontroller;
294 // Phase One: obtain ids: not needed, we have infinite pool
295
296 // Phase Two: add MPQCJobs and send
297 const size_t NoJobs = mpqccontroller.addJobsFromQueue(
298 params.DoLongrange.get() ? MPQCData::DoSampleDensity : MPQCData::DontSampleDensity,
299 params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
300 params.executable.get().string()
301 );
302 if (NoJobs != NumberJobs) {
303 ELOG(1, "Not all jobs were added succesfully.");
304 return Action::failure;
305 }
306 LOG(1, "INFO: Added " << NoJobs << " from FragmentJobsQueue.");
307
308 // prepare process
309 setMaxSteps(NumberJobs);
310 mpqccontroller.setUpdateHandler(
311 boost::bind(&updateSteps, boost::ref(*this), _1, _2)
312 );
313 start();
314 mpqccontroller.run();
315 stop();
316 if (mpqccontroller.getExitflag() != 0) {
317 ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag before getting results.");
318 return Action::failure;
319 }
320
321 // get back the results and place them in shortrangedata
322 mpqccontroller.getResults(shortrangedata);
323 ASSERT( shortrangedata.size() == NumberJobs,
324 "FragmentationFragmentationAutomationAction::performCall() - number of converted results "
325 +toString(shortrangedata.size())+" and number of jobs "+toString(NumberJobs)+ " differ.");
326 if (!checkResults(shortrangedata)) {
327 ELOG(1, "At least one of the fragments' energy could not be properly minimized.");
328 return Action::failure;
329 }
330 if (mpqccontroller.getExitflag() != 0) {
331 ELOG(1, "MPQCCommandFragmentController returned non-zero exit flag after getting results.");
332 return Action::failure;
333 }
334
335 Exitflag += mpqccontroller.getExitflag();
336 }
337#endif
338
339#if defined(HAVE_JOBMARKET) && defined(HAVE_VMG)
340 if (params.DoLongrange.get()) {
341 if ( const_cast<const World &>(World::getInstance()).getAllAtoms().size() == 0) {
342 STATUS("Please load the full molecule into the world before starting this action.");
343 return Action::failure;
344 }
345
346 // obtain combined charge density
347 std::map<JobId_t, MPQCData> shortrangedata_downsampled;
348 SamplingGridProperties domain(ExportGraph_ToJobs::getDomainGrid(params.level.get()));
349 FragmentationChargeDensity summedChargeDensity(
350 shortrangedata);
351 {
352 SamplingGrid zero_globalgrid(domain);
353 summedChargeDensity(
354 shortrangedata,
355 FragmentJobQueue::getInstance().getKeySets(),
356 zero_globalgrid);
357 }
358 const std::vector<SamplingGrid> full_sample = summedChargeDensity.getFullSampledGrid();
359 LOG(1, "INFO: There are " << shortrangedata.size() << " short-range and "
360 << full_sample.size() << " level-wise long-range jobs.");
361
362 // check boundary conditions
363 const BoundaryConditions::Conditions_t &conditions =
364 World::getInstance().getDomain().getConditions();
365 const bool OpenBoundaryConditions =
366 !((conditions[0] == BoundaryConditions::Wrap) &&
367 (conditions[1] == BoundaryConditions::Wrap) &&
368 (conditions[2] == BoundaryConditions::Wrap));
369 LOG(1, std::string("INFO: Using ")
370 << (OpenBoundaryConditions ? "open" : "periodic")
371 << " boundary conditions.");
372
373 // Phase Four: obtain more ids
374 std::map<JobId_t, VMGData> longrangedata;
375 {
376 VMGFragmentController vmgcontroller(io_service);
377 vmgcontroller.setHost(params.host.get());
378 vmgcontroller.setPort(params.port.get());
379 const size_t NoJobs = shortrangedata.size()+full_sample.size();
380 vmgcontroller.requestIds(2*NoJobs);
381 if (vmgcontroller.getExitflag() != 0) {
382 ELOG(1, "VMGFragmentController returned non-zero exit flag on requesting long-range ids.");
383 return Action::failure;
384 }
385
386 // Phase Five a: create VMGJobs for electronic charge distribution
387 const size_t near_field_cells = params.near_field_cells.get();
388 const size_t interpolation_degree = params.interpolation_degree.get();
389 if (!vmgcontroller.createLongRangeJobs(
390 shortrangedata,
391 full_sample,
392 near_field_cells,
393 interpolation_degree,
394 VMGFragmentController::DontSampleParticles,
395 VMGFragmentController::DoTreatGrid,
396 params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
397 params.DoPrintDebug.get(),
398 OpenBoundaryConditions,
399 params.DoSmearCharges.get(),
400 false)) {
401 STATUS("Could not create long-range jobs for electronic charge distribution.");
402 return Action::failure;
403 }
404
405 // Phase Six a: calculate result
406 vmgcontroller.waitforResults(NoJobs);
407 if (vmgcontroller.getExitflag() != 0) {
408 ELOG(1, "VMGFragmentController returned non-zero exit flag before getting results.");
409 return Action::failure;
410 }
411 vmgcontroller.getResults(longrangedata);
412 ASSERT( NoJobs == longrangedata.size(),
413 "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
414 +toString(full_sample.size())+"="+toString(NoJobs)
415 +" and first VMGresults "+toString(longrangedata.size())+" don't match.");
416 if (!checkResults(longrangedata)) {
417 ELOG(1, "At least one of the fragments' electronic long-range potential could not be properly minimized.");
418 return Action::failure;
419 }
420 Exitflag += vmgcontroller.getExitflag();
421
422 {
423 std::map<JobId_t, VMGData> longrangedata_both;
424 // Phase Five b: create VMGJobs for nuclei charge distributions
425 const size_t near_field_cells = params.near_field_cells.get();
426 const size_t interpolation_degree = params.interpolation_degree.get();
427 if (!vmgcontroller.createLongRangeJobs(
428 shortrangedata,
429 full_sample,
430 near_field_cells,
431 interpolation_degree,
432 VMGFragmentController::DoSampleParticles,
433 VMGFragmentController::DoTreatGrid,
434 params.DoValenceOnly.get() ? MPQCData::DoSampleValenceOnly : MPQCData::DontSampleValenceOnly,
435 params.DoPrintDebug.get(),
436 OpenBoundaryConditions,
437 params.DoSmearCharges.get(),
438 params.UseImplicitCharges.get())) {
439 STATUS("Could not create long-range jobs for nuclei charge distribution.");
440 return Action::failure;
441 }
442
443 // Phase Six b: calculate result
444 vmgcontroller.waitforResults(NoJobs);
445 if (vmgcontroller.getExitflag() != 0) {
446 ELOG(1, "VMGFragmentController returned non-zero exit flag before getting results.");
447 return Action::failure;
448 }
449 vmgcontroller.getResults(longrangedata_both);
450 ASSERT( NoJobs == longrangedata_both.size(),
451 "FragmentationFragmentationAutomationAction::performCall() - number of MPQCresults+"
452 +toString(full_sample.size())+"="+toString(NoJobs)
453 +" and second VMGresults "+toString(longrangedata_both.size())+" don't match.");
454 if (!checkResults(longrangedata_both)) {
455 ELOG(1, "At least one of the fragments' nuclei long-range potential could not be properly minimized.");
456 return Action::failure;
457 }
458 if (vmgcontroller.getExitflag() != 0) {
459 ELOG(1, "VMGFragmentController returned non-zero exit flag after getting results.");
460 return Action::failure;
461 }
462 Exitflag += vmgcontroller.getExitflag();
463
464 // now replace nuclei-nuclei contribution with values from direct summation (is exact for
465 // open boundary conditions
466 //
467 // NOTE: We need to calculate both_sampled_potential for fitting partial charges,
468 // nowhere else is this quantity needed.
469 //if (OpenBoundaryConditions)
470 // calculateNucleiNucleiLongRangeContribution(shortrangedata, longrangedata_both);
471
472 // go through either data and replace nuclei_long with contribution from both
473 ASSERT( longrangedata.size() == longrangedata_both.size(),
474 "FragmentationFragmentationAutomationAction::performCall() - longrange results have different sizes.");
475 std::map<JobId_t, VMGData>::iterator destiter = longrangedata.begin();
476 std::map<JobId_t, VMGData>::const_iterator srciter = longrangedata_both.begin();
477 for (;destiter != longrangedata.end(); ++srciter, ++destiter) {
478 destiter->second.both_sampled_potential = srciter->second.sampled_potential;
479 destiter->second.nuclei_long = srciter->second.nuclei_long;
480 destiter->second.forces = srciter->second.forces;
481 destiter->second.hasForces = srciter->second.hasForces;
482 }
483 }
484 }
485
486 if (params.DoPrintDebug.get()) {
487 std::map<JobId_t, std::string> debugData;
488 {
489 if (!full_sample.empty()) {
490 // create debug jobs for each level to print the summed-up potential to vtk files
491 VMGDebugGridFragmentController debugcontroller(io_service);
492 debugcontroller.setHost(params.host.get());
493 debugcontroller.setPort(params.port.get());
494 debugcontroller.requestIds(full_sample.size());
495 if (!debugcontroller.createDebugJobs(full_sample, OpenBoundaryConditions)) {
496 STATUS("Could not create debug jobs.");
497 return Action::failure;
498 }
499 debugcontroller.waitforResults(full_sample.size());
500 debugcontroller.getResults(debugData);
501 Exitflag += debugcontroller.getExitflag();
502 }
503 }
504 }
505 container.clear();
506 container.addFullResults(keysets, forcekeysets, edges_per_fragment, shortrangedata, longrangedata);
507 } else {
508 container.clear();
509 container.addShortRangeResults(keysets, forcekeysets, edges_per_fragment, shortrangedata);
510 }
511#else
512 container.clear();
513 container.addShortRangeResults(keysets, forcekeysets, edges_per_fragment, shortrangedata);
514#endif
515
516 // now clear all present jobs as we are done
517 FragmentJobQueue::getInstance().clear();
518
519 if (Exitflag != 0)
520 STATUS("Controller has returned failure.");
521 return (Exitflag == 0) ? Action::success : Action::failure;
522}
523
524ActionState::ptr FragmentationFragmentationAutomationAction::performUndo(ActionState::ptr _state) {
525 return Action::success;
526}
527
528ActionState::ptr FragmentationFragmentationAutomationAction::performRedo(ActionState::ptr _state){
529 return Action::success;
530}
531
532bool FragmentationFragmentationAutomationAction::canUndo() {
533 return false;
534}
535
536bool FragmentationFragmentationAutomationAction::shouldUndo() {
537 return false;
538}
539/** =========== end of function ====================== */
Note: See TracBrowser for help on using the repository browser.