source: src/LinkedCell/LinkedCell_Model.cpp@ b75386

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

Verbosity fixes to LinkedCell_Model and _View.

  • Property mode set to 100644
File size: 17.7 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2012 University of Bonn. All rights reserved.
5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
8/*
9 * LinkedCell_Model.cpp
10 *
11 * Created on: Nov 15, 2011
12 * Author: heber
13 */
14
15// include config.h
16#ifdef HAVE_CONFIG_H
17#include <config.h>
18#endif
19
20#include "CodePatterns/MemDebug.hpp"
21
22#include "LinkedCell_Model.hpp"
23
24#include <algorithm>
25#include <boost/bind.hpp>
26#include <boost/multi_array.hpp>
27#include <limits>
28
29#include "Atom/AtomObserver.hpp"
30#include "Atom/TesselPoint.hpp"
31#include "Box.hpp"
32#include "CodePatterns/Assert.hpp"
33#include "CodePatterns/Info.hpp"
34#include "CodePatterns/Log.hpp"
35#include "CodePatterns/Observer/Observer.hpp"
36#include "CodePatterns/Observer/Notification.hpp"
37#include "CodePatterns/toString.hpp"
38#include "LinearAlgebra/RealSpaceMatrix.hpp"
39#include "LinearAlgebra/Vector.hpp"
40#include "LinkedCell/IPointCloud.hpp"
41#include "LinkedCell/LinkedCell.hpp"
42#include "LinkedCell/LinkedCell_Model_changeModel.hpp"
43#include "LinkedCell/LinkedCell_Model_LinkedCellArrayCache.hpp"
44#include "World.hpp"
45
46// initialize static entities
47LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::NearestNeighbors;
48
49
50/** Constructor of LinkedCell_Model.
51 *
52 * @param radius desired maximum neighborhood distance
53 * @param _domain Box instance with domain size and boundary conditions
54 */
55LinkedCell::LinkedCell_Model::LinkedCell_Model(const double radius, const Box &_domain) :
56 ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
57 Changes( new changeModel(radius) ),
58 internal_Sizes(NULL),
59 N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
60 domain(_domain)
61{
62 // set default argument
63 NearestNeighbors[0] = NearestNeighbors[1] = NearestNeighbors[2] = 1;
64
65 // get the partition of the domain
66 setPartition(radius);
67
68 // allocate linked cell structure
69 AllocateCells();
70
71 // sign in to AtomObserver
72 startListening();
73}
74
75/** Constructor of LinkedCell_Model.
76 *
77 * @oaram set set of points to place into the linked cell structure
78 * @param radius desired maximum neighborhood distance
79 * @param _domain Box instance with domain size and boundary conditions
80 */
81LinkedCell::LinkedCell_Model::LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain) :
82 ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
83 Changes( new changeModel(radius) ),
84 internal_Sizes(NULL),
85 N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
86 domain(_domain)
87{
88 Info info(__func__);
89
90 // get the partition of the domain
91 setPartition(radius);
92
93 // allocate linked cell structure
94 AllocateCells();
95
96 insertPointCloud(set);
97
98 // sign in to AtomObserver
99 startListening();
100}
101
102/** Destructor of class LinkedCell_Model.
103 *
104 */
105LinkedCell::LinkedCell_Model::~LinkedCell_Model()
106{
107 // sign off from observables
108 stopListening();
109
110 // delete change queue
111 delete Changes;
112
113 // reset linked cell structure
114 Reset();
115
116}
117
118/** Signs in to AtomObserver and World to known about all changes.
119 *
120 */
121void LinkedCell::LinkedCell_Model::startListening()
122{
123 World::getInstance().signOn(this, World::AtomInserted);
124 World::getInstance().signOn(this, World::AtomRemoved);
125 AtomObserver::getInstance().signOn(this, AtomObservable::PositionChanged);
126}
127
128/** Signs off from AtomObserver and World.
129 *
130 */
131void LinkedCell::LinkedCell_Model::stopListening()
132{
133 World::getInstance().signOff(this, World::AtomInserted);
134 World::getInstance().signOff(this, World::AtomRemoved);
135 AtomObserver::getInstance().signOff(this, AtomObservable::PositionChanged);
136}
137
138
139/** Allocates as much cells per axis as required by
140 * LinkedCell_Model::BoxPartition.
141 *
142 */
143void LinkedCell::LinkedCell_Model::AllocateCells()
144{
145 // resize array
146 tripleIndex index;
147 for (int i=0;i<NDIM;i++)
148 index[i] = static_cast<LinkedCellArray::index>(Dimensions.at(i,i));
149 N->setN().resize(index);
150 ASSERT(getSize(0)*getSize(1)*getSize(2) < MAX_LINKEDCELLNODES,
151 "LinkedCell_Model::AllocateCells() - Number linked of linked cell nodes exceeded hard-coded limit, use greater edge length!");
152 LOG(1, "INFO: Allocating array ("
153 << getSize(0) << ","
154 << getSize(1) << ","
155 << getSize(2) << ") for a new LinkedCell_Model.");
156
157 // allocate LinkedCell instances
158 for(index[0] = 0; index[0] != static_cast<LinkedCellArray::index>(Dimensions.at(0,0)); ++index[0]) {
159 for(index[1] = 0; index[1] != static_cast<LinkedCellArray::index>(Dimensions.at(1,1)); ++index[1]) {
160 for(index[2] = 0; index[2] != static_cast<LinkedCellArray::index>(Dimensions.at(2,2)); ++index[2]) {
161 LOG(5, "INFO: Creating cell at " << index[0] << " " << index[1] << " " << index[2] << ".");
162 N->setN()(index) = new LinkedCell(index);
163 }
164 }
165 }
166}
167
168/** Frees all Linked Cell instances and sets array dimensions to (0,0,0).
169 *
170 */
171void LinkedCell::LinkedCell_Model::Reset()
172{
173 // free all LinkedCell instances
174 for(iterator3 iter3 = N->setN().begin(); iter3 != N->setN().end(); ++iter3) {
175 for(iterator2 iter2 = (*iter3).begin(); iter2 != (*iter3).end(); ++iter2) {
176 for(iterator1 iter1 = (*iter2).begin(); iter1 != (*iter2).end(); ++iter1) {
177 delete *iter1;
178 }
179 }
180 }
181 // set dimensions to zero
182 N->setN().resize(boost::extents[0][0][0]);
183}
184
185/** Inserts all points contained in \a set.
186 *
187 * @param set set with points to insert into linked cell structure
188 */
189void LinkedCell::LinkedCell_Model::insertPointCloud(IPointCloud &set)
190{
191 if (set.IsEmpty()) {
192 ELOG(1, "set is NULL or contains no linked cell nodes!");
193 return;
194 }
195
196 // put each atom into its respective cell
197 set.GoToFirst();
198 while (!set.IsEnd()) {
199 TesselPoint *Walker = set.GetPoint();
200 addNode(Walker);
201 set.GoToNext();
202 }
203}
204
205/** Calculates the required edge length for the given desired distance.
206 *
207 * We need to make some matrix transformations in order to obtain the required
208 * edge lengths per axis. Goal is guarantee that whatever the shape of the
209 * domain that always all points at least up to \a distance away are contained
210 * in the nearest neighboring cells.
211 *
212 * @param distance distance of this linked cell array
213 */
214void LinkedCell::LinkedCell_Model::setPartition(double distance)
215{
216 // generate box matrix of desired edge length
217 RealSpaceMatrix neighborhood;
218 neighborhood.setIdentity();
219 neighborhood *= distance;
220
221 // obtain refs to both domain matrix transformations
222 //const RealSpaceMatrix &M = domain.getM();
223 const RealSpaceMatrix &Minv = domain.getMinv();
224
225 RealSpaceMatrix output = Minv * neighborhood;
226
227 //std::cout << Minv << " * " << neighborhood << " = " << output << std::endl;
228
229 Dimensions = output.invert();
230 Partition = Minv*Dimensions; //
231
232 //std::cout << "Dimensions are then " << Dimensions << std::endl;
233 //std::cout << "Partition matrix is then " << Partition << std::endl;
234}
235
236/** Returns the number of required neighbor-shells to get all neighboring points
237 * in the given \a distance.
238 *
239 * @param distance radius of neighborhood sphere
240 * @return number of LinkedCell's per dimension to get all neighbors
241 */
242const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getStep(const double distance) const
243{
244 tripleIndex index;
245 index[0] = index[1] = index[2] = 0;
246
247 if (fabs(distance) < std::numeric_limits<double>::min())
248 return index;
249 // generate box matrix of desired edge length
250 RealSpaceMatrix neighborhood;
251 neighborhood.setIdentity();
252 neighborhood *= distance;
253
254 const RealSpaceMatrix output = Partition * neighborhood;
255
256 //std::cout << "GetSteps: " << Partition << " * " << neighborhood << " = " << output << std::endl;
257
258 const RealSpaceMatrix steps = output;
259 for (size_t i =0; i<NDIM; ++i)
260 index[i] = ceil(steps.at(i,i));
261 LOG(2, "INFO: number of shells are ("+toString(index[0])+","+toString(index[1])+","+toString(index[2])+").");
262
263 return index;
264}
265
266/** Calculates the index of the cell \a position would belong to.
267 *
268 * @param position position whose associated cell to calculate
269 * @return index of the cell
270 */
271const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
272{
273 tripleIndex index;
274 Vector x(Partition*position);
275 LOG(2, "INFO: Transformed position is " << x << ".");
276 for (int i=0;i<NDIM;i++) {
277 index[i] = static_cast<LinkedCellArray::index>(floor(x[i]));
278 }
279 return index;
280}
281
282/** Adds an update to the list of lazy changes to add a node.
283 *
284 * @param Walker node to add
285 */
286void LinkedCell::LinkedCell_Model::addNode(const TesselPoint *Walker)
287{
288 LOG(2, "INFO: Requesting update to add node " << *Walker << ".");
289 Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::addNode_internal, this, _1));
290}
291
292/** Adds an update to the list of lazy changes to add remove a node.
293 *
294 * We do nothing of Walker is not found
295 *
296 * @param Walker node to remove
297 */
298void LinkedCell::LinkedCell_Model::deleteNode(const TesselPoint *Walker)
299{
300 LOG(2, "INFO: Requesting update to delete node " << *Walker << ".");
301 Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::deleteNode_internal, this, _1));
302}
303
304/** Adds an update to the list of lazy changes to move a node.
305 *
306 * @param Walker node who has moved.
307 */
308void LinkedCell::LinkedCell_Model::moveNode(const TesselPoint *Walker)
309{
310 LOG(2, "INFO: Requesting update to move node " << *Walker << " to position "
311 << Walker->getPosition() << ".");
312 Changes->addUpdate(Walker, 10, boost::bind(&LinkedCell_Model::moveNode_internal, this, _1));
313}
314
315/** Internal function to add a node to the linked cell structure
316 *
317 * @param Walker node to add
318 */
319void LinkedCell::LinkedCell_Model::addNode_internal(const TesselPoint *Walker)
320{
321 // find index
322 tripleIndex index = getIndexToVector(Walker->getPosition());
323 LOG(2, "INFO: " << *Walker << " goes into cell " << index[0] << ", " << index[1] << ", " << index[2] << ".");
324 LOG(2, "INFO: Cell's indices are "
325 << (N->getN())(index)->getIndex(0) << " "
326 << (N->getN())(index)->getIndex(1) << " "
327 << (N->getN())(index)->getIndex(2) << ".");
328 // add to cell
329 (N->setN())(index)->addPoint(Walker);
330 // add to index with check for presence
331 std::pair<MapPointToCell::iterator, bool> inserter = CellLookup.insert( std::make_pair(Walker, (N->setN())(index)) );
332 ASSERT( inserter.second,
333 "LinkedCell_Model::addNode() - Walker "
334 +toString(*Walker)+" is already present in cell "
335 +toString((inserter.first)->second->getIndex(0))+" "
336 +toString((inserter.first)->second->getIndex(1))+" "
337 +toString((inserter.first)->second->getIndex(2))+".");
338}
339
340/** Internal function to remove a node to the linked cell structure
341 *
342 * We do nothing of Walker is not found
343 *
344 * @param Walker node to remove
345 */
346void LinkedCell::LinkedCell_Model::deleteNode_internal(const TesselPoint *Walker)
347{
348 MapPointToCell::iterator iter = CellLookup.find(Walker);
349 ASSERT(iter != CellLookup.end(),
350 "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
351 if (iter != CellLookup.end()) {
352 // remove from lookup
353 CellLookup.erase(iter);
354 // remove from cell
355 iter->second->deletePoint(Walker);
356 }
357}
358
359/** Internal function to move node from current cell to another on position change.
360 *
361 * @param Walker node who has moved.
362 */
363void LinkedCell::LinkedCell_Model::moveNode_internal(const TesselPoint *Walker)
364{
365 MapPointToCell::iterator iter = CellLookup.find(Walker);
366 ASSERT(iter != CellLookup.end(),
367 "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
368 if (iter != CellLookup.end()) {
369 tripleIndex index = getIndexToVector(Walker->getPosition());
370 if (index != iter->second->getIndices()) {
371 // remove in old cell
372 iter->second->deletePoint(Walker);
373 // add to new cell
374 N->setN()(index)->addPoint(Walker);
375 // update lookup
376 iter->second = N->setN()(index);
377 }
378 }
379}
380
381/** Checks whether cell indicated by \a relative relative to LinkedCell_Model::internal_index
382 * is out of bounds.
383 *
384 * \note We do not check for boundary conditions of LinkedCell_Model::domain,
385 * we only look at the array sizes
386 *
387 * @param relative index relative to LinkedCell_Model::internal_index.
388 * @return true - relative index is still inside bounds, false - outside
389 */
390bool LinkedCell::LinkedCell_Model::checkArrayBounds(const tripleIndex &index) const
391{
392 bool status = true;
393 for (size_t i=0;i<NDIM;++i) {
394 status = status && (
395 (index[i] >= 0) &&
396 (index[i] < getSize(i))
397 );
398 }
399 return status;
400}
401
402/** Corrects \a index according to boundary conditions of LinkedCell_Model::domain .
403 *
404 * @param index index to correct according to boundary conditions
405 */
406void LinkedCell::LinkedCell_Model::applyBoundaryConditions(tripleIndex &index) const
407{
408 for (size_t i=0;i<NDIM;++i) {
409 switch (domain.getConditions()[i]) {
410 case Box::Wrap:
411 if ((index[i] < 0) || (index[i] >= getSize(i)))
412 index[i] = (index[i] % getSize(i));
413 break;
414 case Box::Bounce:
415 if (index[i] < 0)
416 index[i] = 0;
417 if (index[i] >= getSize(i))
418 index[i] = getSize(i)-1;
419 break;
420 case Box::Ignore:
421 if (index[i] < 0)
422 index[i] = 0;
423 if (index[i] >= getSize(i))
424 index[i] = getSize(i)-1;
425 break;
426 default:
427 ASSERT(0, "LinkedCell_Model::checkBounds() - unknown boundary conditions.");
428 break;
429 }
430 }
431}
432
433/** Calculates the interval bounds of the linked cell grid.
434 *
435 * \note we assume for index to allows be valid, i.e. within the range of LinkedCell_Model::N.
436 *
437 * \param index index to give relative bounds to
438 * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
439 * \return pair of tripleIndex indicating lower and upper bounds
440 */
441const LinkedCell::LinkedCell_Model::LinkedCellNeighborhoodBounds LinkedCell::LinkedCell_Model::getNeighborhoodBounds(
442 const tripleIndex &index,
443 const tripleIndex &step
444 ) const
445{
446 LinkedCellNeighborhoodBounds neighbors;
447
448 // calculate bounds
449 for (size_t i=0;i<NDIM;++i) {
450 ASSERT(index[i] >= 0,
451 "LinkedCell_Model::getNeighborhoodBounds() - index "+toString(index)+" out of lower bounds.");
452 ASSERT (index[i] < getSize(i),
453 "LinkedCell_Model::getNeighborhoodBounds() - index "+toString(index)+" out of upper bounds.");
454 switch (domain.getConditions()[i]) {
455 case Box::Wrap:
456 if ((index[i] - step[i]) < 0)
457 neighbors.first[i] = getSize(i) + (index[i] - step[i]);
458 else if ((index[i] - step[i]) >= getSize(i))
459 neighbors.first[i] = (index[i] - step[i]) - getSize(i);
460 else
461 neighbors.first[i] = index[i] - step[i];
462 neighbors.second[i] = 2*step[i]+1;
463 break;
464 case Box::Bounce:
465 neighbors.second[i] = 2*step[i]+1;
466 if (index[i] - step[i] >= 0) {
467 neighbors.first[i] = index[i] - step[i];
468 } else {
469 neighbors.first[i] = 0;
470 neighbors.second[i] = index[i] + step[i]+1;
471 }
472 if (index[i] + step[i] >= getSize(i)) {
473 neighbors.second[i] = getSize(i) - (index[i] - step[i]);
474 }
475 break;
476 case Box::Ignore:
477 if (index[i] - step[i] < 0)
478 neighbors.first[i] = 0;
479 else
480 neighbors.first[i] = index[i] - step[i];
481 if (index[i] + step[i] >= getSize(i))
482 neighbors.second[i] = getSize(i) - index[i];
483 else
484 neighbors.second[i] = 2*step[i]+1;
485 break;
486 default:
487 ASSERT(0, "LinkedCell_Model::getNeighborhoodBounds() - unknown boundary conditions.");
488 break;
489 }
490 }
491
492 return neighbors;
493}
494
495/** Returns a reference to the cell indicated by LinkedCell_Model::internal_index.
496 *
497 * \return LinkedCell ref to current cell
498 */
499const LinkedCell::LinkedCell& LinkedCell::LinkedCell_Model::getCell(const tripleIndex &index) const
500{
501 return *(N->getN()(index));
502}
503
504
505/** Returns size of array for given \a dim.
506 *
507 * @param dim desired dimension
508 * @return size of array along dimension
509 */
510LinkedCell::LinkedCellArray::index LinkedCell::LinkedCell_Model::getSize(const size_t dim) const
511{
512 ASSERT((dim >= 0) && (dim < NDIM),
513 "LinkedCell_Model::getSize() - dimension "
514 +toString(dim)+" is out of bounds.");
515 return N->getN().shape()[dim];
516}
517
518/** Callback function for Observer mechanism.
519 *
520 * @param publisher reference to the Observable that calls
521 */
522void LinkedCell::LinkedCell_Model::update(Observable *publisher)
523{
524 ELOG(2, "LinkedCell_Model received inconclusive general update from "
525 << publisher << ".");
526}
527
528/** Callback function for the Notifications mechanism.
529 *
530 * @param publisher reference to the Observable that calls
531 * @param notification specific notification as cause of the call
532 */
533void LinkedCell::LinkedCell_Model::recieveNotification(Observable *publisher, Notification_ptr notification)
534{
535 // either it's the World or from the atom (through relay) itself
536 if (publisher == World::getPointer()) {
537 switch(notification->getChannelNo()) {
538 case World::AtomInserted:
539 addNode(World::getInstance().lastChanged<atom>());
540 break;
541 case World::AtomRemoved:
542 deleteNode(World::getInstance().lastChanged<atom>());
543 break;
544 }
545 } else {
546 switch(notification->getChannelNo()) {
547 case AtomObservable::PositionChanged:
548 {
549 moveNode(dynamic_cast<const TesselPoint *>(publisher));
550 break;
551 }
552 default:
553 LOG(2, "LinkedCell_Model received unwanted notification from AtomObserver's channel "
554 << notification->getChannelNo() << ".");
555 break;
556 }
557 }
558}
559
560/** Callback function when an Observer dies.
561 *
562 * @param publisher reference to the Observable that calls
563 */
564void LinkedCell::LinkedCell_Model::subjectKilled(Observable *publisher)
565{}
566
Note: See TracBrowser for help on using the repository browser.