source: src/Box.cpp@ dd067a

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

Updated all source files's copyright note to current year 2012.

  • Property mode set to 100644
File size: 10.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 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
8/*
9 * Box.cpp
10 *
11 * Created on: Jun 30, 2010
12 * Author: crueger
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 "Box.hpp"
23
24#include <cmath>
25#include <iostream>
26#include <cstdlib>
27
28#include "CodePatterns/Assert.hpp"
29#include "CodePatterns/Log.hpp"
30#include "CodePatterns/Observer/Channels.hpp"
31#include "CodePatterns/Observer/Notification.hpp"
32#include "CodePatterns/Verbose.hpp"
33#include "Helpers/defs.hpp"
34#include "LinearAlgebra/RealSpaceMatrix.hpp"
35#include "LinearAlgebra/Vector.hpp"
36#include "LinearAlgebra/Plane.hpp"
37#include "Shapes/BaseShapes.hpp"
38#include "Shapes/ShapeOps.hpp"
39
40
41Box::Box() :
42 Observable("Box"),
43 M(new RealSpaceMatrix()),
44 Minv(new RealSpaceMatrix())
45{
46 internal_list.reserve(pow(3,3));
47 coords.reserve(NDIM);
48 index.reserve(NDIM);
49
50 // observable stuff
51 Channels *OurChannel = new Channels;
52 NotificationChannels.insert( std::make_pair(this, OurChannel) );
53 // add instance for each notification type
54 for (size_t type = 0; type < NotificationType_MAX; ++type)
55 OurChannel->addChannel(type);
56
57 M->setIdentity();
58 Minv->setIdentity();
59 conditions.resize(3);
60 conditions[0] = conditions[1] = conditions[2] = Wrap;
61}
62
63Box::Box(const Box& src) :
64 Observable("Box"),
65 conditions(src.conditions),
66 M(new RealSpaceMatrix(*src.M)),
67 Minv(new RealSpaceMatrix(*src.Minv))
68{
69 internal_list.reserve(pow(3,3));
70 coords.reserve(NDIM);
71 index.reserve(NDIM);
72
73 // observable stuff
74 Channels *OurChannel = new Channels;
75 NotificationChannels.insert( std::make_pair(this, OurChannel) );
76 // add instance for each notification type
77 for (size_t type = 0; type < NotificationType_MAX; ++type)
78 OurChannel->addChannel(type);
79}
80
81Box::Box(RealSpaceMatrix _M) :
82 Observable("Box"),
83 M(new RealSpaceMatrix(_M)),
84 Minv(new RealSpaceMatrix())
85{
86 internal_list.reserve(pow(3,3));
87 coords.reserve(NDIM);
88 index.reserve(NDIM);
89
90 // observable stuff
91 Channels *OurChannel = new Channels;
92 NotificationChannels.insert( std::make_pair(this, OurChannel) );
93 // add instance for each notification type
94 for (size_t type = 0; type < NotificationType_MAX; ++type)
95 OurChannel->addChannel(type);
96
97 ASSERT(M->determinant()!=0,"Matrix in Box construction was not invertible");
98 *Minv = M->invert();
99 conditions.resize(3);
100 conditions[0] = conditions[1] = conditions[2] = Wrap;
101}
102
103Box::~Box()
104{
105 // observable stuff
106 std::map<Observable *, Channels*>::iterator iter = NotificationChannels.find(this);
107 delete iter->second;
108 NotificationChannels.erase(iter);
109
110 delete M;
111 delete Minv;
112}
113
114const RealSpaceMatrix &Box::getM() const{
115 return *M;
116}
117const RealSpaceMatrix &Box::getMinv() const{
118 return *Minv;
119}
120
121void Box::setM(RealSpaceMatrix _M){
122 ASSERT(_M.determinant()!=0,"Matrix in Box construction was not invertible");
123 OBSERVE;
124 NOTIFY(MatrixChanged);
125 *M =_M;
126 *Minv = M->invert();
127}
128
129Vector Box::translateIn(const Vector &point) const{
130 return (*M) * point;
131}
132
133Vector Box::translateOut(const Vector &point) const{
134 return (*Minv) * point;
135}
136
137Vector Box::WrapPeriodically(const Vector &point) const{
138 Vector helper = translateOut(point);
139 for(int i=NDIM;i--;){
140
141 switch(conditions[i]){
142 case Wrap:
143 helper.at(i)=fmod(helper.at(i),1);
144 helper.at(i)+=(helper.at(i)>=0)?0:1;
145 break;
146 case Bounce:
147 {
148 // there probably is a better way to handle this...
149 // all the fabs and fmod modf probably makes it very slow
150 double intpart,fracpart;
151 fracpart = modf(fabs(helper.at(i)),&intpart);
152 helper.at(i) = fabs(fracpart-fmod(intpart,2));
153 }
154 break;
155 case Ignore:
156 break;
157 default:
158 ASSERT(0,"No default case for this");
159 break;
160 }
161
162 }
163 return translateIn(helper);
164}
165
166bool Box::isInside(const Vector &point) const
167{
168 bool result = true;
169 Vector tester = translateOut(point);
170
171 for(int i=0;i<NDIM;i++)
172 result = result &&
173 ((conditions[i] == Ignore) ||
174 ((tester[i] >= -MYEPSILON) &&
175 ((tester[i] - 1.) < MYEPSILON)));
176
177 return result;
178}
179
180
181VECTORSET(std::vector) Box::explode(const Vector &point,int n) const{
182 ASSERT(isInside(point),"Exploded point not inside Box");
183 internal_explode(point, n);
184 VECTORSET(std::vector) res(internal_list);
185 return res;
186}
187
188void Box::internal_explode(const Vector &point,int n) const{
189 internal_list.clear();
190 size_t list_index = 0;
191
192 Vector translater = translateOut(point);
193 Vector mask; // contains the ignored coordinates
194
195 // count the number of coordinates we need to do
196 int dims = 0; // number of dimensions that are not ignored
197 coords.clear();
198 index.clear();
199 for(int i=0;i<NDIM;++i){
200 if(conditions[i]==Ignore){
201 mask[i]=translater[i];
202 continue;
203 }
204 coords.push_back(i);
205 index.push_back(-n);
206 dims++;
207 } // there are max vectors in total we need to create
208 internal_list.resize(pow(2*n+1, dims));
209
210 if(!dims){
211 // all boundaries are ignored
212 internal_list[list_index++] = point;
213 return;
214 }
215
216 bool done = false;
217 while(!done){
218 // create this vector
219 Vector helper;
220 for(int i=0;i<dims;++i){
221 switch(conditions[coords[i]]){
222 case Wrap:
223 helper[coords[i]] = index[i]+translater[coords[i]];
224 break;
225 case Bounce:
226 {
227 // Bouncing the coordinate x produces the series:
228 // 0 -> x
229 // 1 -> 2-x
230 // 2 -> 2+x
231 // 3 -> 4-x
232 // 4 -> 4+x
233 // the first number is the next bigger even number (n+n%2)
234 // the next number is the value with alternating sign (x-2*(n%2)*x)
235 // the negative numbers produce the same sequence reversed and shifted
236 int n = abs(index[i]) + ((index[i]<0)?-1:0);
237 int sign = (index[i]<0)?-1:+1;
238 int even = n%2;
239 helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
240 helper[coords[i]]*=sign;
241 }
242 break;
243 case Ignore:
244 ASSERT(0,"Ignored coordinate handled in generation loop");
245 break;
246 default:
247 ASSERT(0,"No default case for this switch-case");
248 break;
249 }
250
251 }
252 // add back all ignored coordinates (not handled in above loop)
253 helper+=mask;
254 ASSERT(list_index < internal_list.size(),
255 "Box::internal_explode() - we have estimated the number of vectors wrong: "
256 +toString(list_index) +" >= "+toString(internal_list.size())+".");
257 internal_list[list_index++] = translateIn(helper);
258 // set the new indexes
259 int pos=0;
260 ++index[pos];
261 while(index[pos]>n){
262 index[pos++]=-n;
263 if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
264 done = true;
265 break;
266 }
267 ++index[pos];
268 }
269 }
270}
271
272VECTORSET(std::vector) Box::explode(const Vector &point) const{
273 ASSERT(isInside(point),"Exploded point not inside Box");
274 return explode(point,1);
275}
276
277double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
278 Vector helper1(!isInside(point1) ? WrapPeriodically(point1) : point1);
279 Vector helper2(!isInside(point2) ? WrapPeriodically(point2) : point2);
280 internal_explode(helper1,1);
281 double res = internal_list.minDistSquared(helper2);
282 return res;
283}
284
285double Box::periodicDistance(const Vector &point1,const Vector &point2) const{
286 double res;
287 res = sqrt(periodicDistanceSquared(point1,point2));
288 return res;
289}
290
291double Box::DistanceToBoundary(const Vector &point) const
292{
293 std::map<double, Plane> DistanceSet;
294 std::vector<std::pair<Plane,Plane> > Boundaries = getBoundingPlanes();
295 for (int i=0;i<NDIM;++i) {
296 const double tempres1 = Boundaries[i].first.distance(point);
297 const double tempres2 = Boundaries[i].second.distance(point);
298 DistanceSet.insert( make_pair(tempres1, Boundaries[i].first) );
299 LOG(1, "Inserting distance " << tempres1 << " and " << tempres2 << ".");
300 DistanceSet.insert( make_pair(tempres2, Boundaries[i].second) );
301 }
302 ASSERT(!DistanceSet.empty(), "Box::DistanceToBoundary() - no distances in map!");
303 return (DistanceSet.begin())->first;
304}
305
306Shape Box::getShape() const{
307 return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
308}
309
310const Box::Conditions_t Box::getConditions() const
311{
312 return conditions;
313}
314
315void Box::setCondition(int i,Box::BoundaryCondition_t condition)
316{
317 OBSERVE;
318 NOTIFY(BoundaryConditionsChanged);
319 conditions[i]=condition;
320}
321
322const std::vector<std::pair<Plane,Plane> > Box::getBoundingPlanes() const
323{
324 std::vector<std::pair<Plane,Plane> > res;
325 for(int i=0;i<NDIM;++i){
326 Vector base1,base2,base3;
327 base2[(i+1)%NDIM] = 1.;
328 base3[(i+2)%NDIM] = 1.;
329 Plane p1(translateIn(base1),
330 translateIn(base2),
331 translateIn(base3));
332 Vector offset;
333 offset[i]=1;
334 Plane p2(translateIn(base1+offset),
335 translateIn(base2+offset),
336 translateIn(base3+offset));
337 res.push_back(make_pair(p1,p2));
338 }
339 ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!");
340 return res;
341}
342
343void Box::setCuboid(const Vector &endpoint)
344{
345 OBSERVE;
346 NOTIFY(MatrixChanged);
347 ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
348 M->setIdentity();
349 M->diagonal()=endpoint;
350 Vector &dinv = Minv->diagonal();
351 for(int i=NDIM;i--;)
352 dinv[i]=1/endpoint[i];
353}
354
355Box &Box::operator=(const Box &src)
356{
357 if(&src!=this){
358 OBSERVE;
359 // new matrix
360 NOTIFY(MatrixChanged);
361 delete M;
362 delete Minv;
363 M = new RealSpaceMatrix(*src.M);
364 Minv = new RealSpaceMatrix(*src.Minv);
365 // new boundary conditions
366 NOTIFY(BoundaryConditionsChanged);
367 conditions = src.conditions;
368 }
369 return *this;
370}
371
372Box &Box::operator=(const RealSpaceMatrix &mat)
373{
374 OBSERVE;
375 NOTIFY(MatrixChanged);
376 setM(mat);
377 return *this;
378}
379
380std::ostream & operator << (std::ostream& ost, const Box &m)
381{
382 ost << m.getM();
383 return ost;
384}
Note: See TracBrowser for help on using the repository browser.