source: src/Box.cpp@ fae462

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 fae462 was 25aa5a, checked in by Frederik Heber <heber@…>, 10 years ago

FIX: Box::setM() did not use const ref as argument.

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