source: src/boundary.cpp@ cd7a87

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 cd7a87 was 9df680, checked in by Frederik Heber <heber@…>, 14 years ago

FIX: FillVoidWithMolecule() now centers filler molecule at its center of gravity.

  • this is the only sensible option if we later rotate the molecule.
  • atom::clone() reduced a lot, as copy constructors contains it all.
  • FillVoidWithMolecule() - filler is now selected to reside entirely within the domain. Otherwise, we seg'fault on save on exit, as one of the father atoms for all copied atoms is missing and hence, we cannot look at its additional...Data for various parsers.
  • new function molecule::removeAtomsInMolecule() which destroys all the molecule's atoms.

TESTFIXES:

  • Filling/3: filled water box has changed due to different center of water molecules.
  • Property mode set to 100644
File size: 61.3 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010 University of Bonn. All rights reserved.
5 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
6 */
7
8/** \file boundary.cpp
9 *
10 * Implementations and super-function for envelopes
11 */
12
13// include config.h
14#ifdef HAVE_CONFIG_H
15#include <config.h>
16#endif
17
18#include "CodePatterns/MemDebug.hpp"
19
20#include "Actions/MoleculeAction/RotateToPrincipalAxisSystemAction.hpp"
21#include "BoundaryPointSet.hpp"
22#include "BoundaryLineSet.hpp"
23#include "BoundaryTriangleSet.hpp"
24#include "CandidateForTesselation.hpp"
25//#include "TesselPoint.hpp"
26#include "World.hpp"
27#include "atom.hpp"
28#include "bond.hpp"
29#include "boundary.hpp"
30#include "config.hpp"
31#include "element.hpp"
32#include "Helpers/helpers.hpp"
33#include "CodePatterns/Info.hpp"
34#include "linkedcell.hpp"
35#include "CodePatterns/Verbose.hpp"
36#include "CodePatterns/Log.hpp"
37#include "molecule.hpp"
38#include "tesselation.hpp"
39#include "tesselationhelpers.hpp"
40#include "World.hpp"
41#include "LinearAlgebra/Plane.hpp"
42#include "LinearAlgebra/RealSpaceMatrix.hpp"
43#include "Box.hpp"
44
45#include <iostream>
46#include <iomanip>
47
48#include<gsl/gsl_poly.h>
49#include<time.h>
50
51// ========================================== F U N C T I O N S =================================
52
53
54/** Determines greatest diameters of a cluster defined by its convex envelope.
55 * Looks at lines parallel to one axis and where they intersect on the projected planes
56 * \param *out output stream for debugging
57 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
58 * \param *mol molecule structure representing the cluster
59 * \param *&TesselStruct Tesselation structure with triangles
60 * \param IsAngstroem whether we have angstroem or atomic units
61 * \return NDIM array of the diameters
62 */
63double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
64{
65 Info FunctionInfo(__func__);
66 // get points on boundary of NULL was given as parameter
67 bool BoundaryFreeFlag = false;
68 double OldComponent = 0.;
69 double tmp = 0.;
70 double w1 = 0.;
71 double w2 = 0.;
72 Vector DistanceVector;
73 Vector OtherVector;
74 int component = 0;
75 int Othercomponent = 0;
76 Boundaries::const_iterator Neighbour;
77 Boundaries::const_iterator OtherNeighbour;
78 double *GreatestDiameter = new double[NDIM];
79
80 const Boundaries *BoundaryPoints;
81 if (BoundaryPtr == NULL) {
82 BoundaryFreeFlag = true;
83 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
84 } else {
85 BoundaryPoints = BoundaryPtr;
86 DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
87 }
88 // determine biggest "diameter" of cluster for each axis
89 for (int i = 0; i < NDIM; i++)
90 GreatestDiameter[i] = 0.;
91 for (int axis = 0; axis < NDIM; axis++)
92 { // regard each projected plane
93 //Log() << Verbose(1) << "Current axis is " << axis << "." << endl;
94 for (int j = 0; j < 2; j++)
95 { // and for both axis on the current plane
96 component = (axis + j + 1) % NDIM;
97 Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
98 //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
99 for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
100 //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl;
101 // seek for the neighbours pair where the Othercomponent sign flips
102 Neighbour = runner;
103 Neighbour++;
104 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
105 Neighbour = BoundaryPoints[axis].begin();
106 DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
107 do { // seek for neighbour pair where it flips
108 OldComponent = DistanceVector[Othercomponent];
109 Neighbour++;
110 if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
111 Neighbour = BoundaryPoints[axis].begin();
112 DistanceVector = (runner->second.second->getPosition()) - (Neighbour->second.second->getPosition());
113 //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
114 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
115 OldComponent) - DistanceVector[Othercomponent] / fabs(
116 DistanceVector[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
117 if (runner != Neighbour) {
118 OtherNeighbour = Neighbour;
119 if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
120 OtherNeighbour = BoundaryPoints[axis].end();
121 OtherNeighbour--;
122 //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
123 // now we have found the pair: Neighbour and OtherNeighbour
124 OtherVector = (runner->second.second->getPosition()) - (OtherNeighbour->second.second->getPosition());
125 //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
126 //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
127 // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
128 w1 = fabs(OtherVector[Othercomponent]);
129 w2 = fabs(DistanceVector[Othercomponent]);
130 tmp = fabs((w1 * DistanceVector[component] + w2
131 * OtherVector[component]) / (w1 + w2));
132 // mark if it has greater diameter
133 //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
134 GreatestDiameter[component] = (GreatestDiameter[component]
135 > tmp) ? GreatestDiameter[component] : tmp;
136 } //else
137 //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
138 }
139 }
140 }
141 Log() << Verbose(0) << "RESULT: The biggest diameters are "
142 << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
143 << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
144 : "atomiclength") << "." << endl;
145
146 // free reference lists
147 if (BoundaryFreeFlag)
148 delete[] (BoundaryPoints);
149
150 return GreatestDiameter;
151}
152;
153
154
155/** Determines the boundary points of a cluster.
156 * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
157 * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
158 * center and first and last point in the triple, it is thrown out.
159 * \param *out output stream for debugging
160 * \param *mol molecule structure representing the cluster
161 * \param *&TesselStruct pointer to Tesselation structure
162 */
163Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
164{
165 Info FunctionInfo(__func__);
166 PointMap PointsOnBoundary;
167 LineMap LinesOnBoundary;
168 TriangleMap TrianglesOnBoundary;
169 Vector *MolCenter = mol->DetermineCenterOfAll();
170 Vector helper;
171 BoundariesTestPair BoundaryTestPair;
172 Vector AxisVector;
173 Vector AngleReferenceVector;
174 Vector AngleReferenceNormalVector;
175 Vector ProjectedVector;
176 Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
177 double angle = 0.;
178
179 // 3a. Go through every axis
180 for (int axis = 0; axis < NDIM; axis++) {
181 AxisVector.Zero();
182 AngleReferenceVector.Zero();
183 AngleReferenceNormalVector.Zero();
184 AxisVector[axis] = 1.;
185 AngleReferenceVector[(axis + 1) % NDIM] = 1.;
186 AngleReferenceNormalVector[(axis + 2) % NDIM] = 1.;
187
188 DoLog(1) && (Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl);
189
190 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
191 for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
192 ProjectedVector = (*iter)->getPosition() - (*MolCenter);
193 ProjectedVector.ProjectOntoPlane(AxisVector);
194
195 // correct for negative side
196 const double radius = ProjectedVector.NormSquared();
197 if (fabs(radius) > MYEPSILON)
198 angle = ProjectedVector.Angle(AngleReferenceVector);
199 else
200 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
201
202 //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
203 if (ProjectedVector.ScalarProduct(AngleReferenceNormalVector) > 0) {
204 angle = 2. * M_PI - angle;
205 }
206 DoLog(1) && (Log() << Verbose(1) << "Inserting " << **iter << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
207 BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, TesselPointDistancePair (radius, (*iter))));
208 if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
209 DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl);
210 DoLog(2) && (Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl);
211 DoLog(2) && (Log() << Verbose(2) << "New vector: " << **iter << endl);
212 const double ProjectedVectorNorm = ProjectedVector.NormSquared();
213 if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
214 BoundaryTestPair.first->second.first = ProjectedVectorNorm;
215 BoundaryTestPair.first->second.second = (*iter);
216 DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
217 } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
218 helper = (*iter)->getPosition() - (*MolCenter);
219 const double oldhelperNorm = helper.NormSquared();
220 helper = BoundaryTestPair.first->second.second->getPosition() - (*MolCenter);
221 if (helper.NormSquared() < oldhelperNorm) {
222 BoundaryTestPair.first->second.second = (*iter);
223 DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
224 } else {
225 DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl);
226 }
227 } else {
228 DoLog(2) && (Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
229 }
230 }
231 }
232 // printing all inserted for debugging
233 // {
234 // Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
235 // int i=0;
236 // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
237 // if (runner != BoundaryPoints[axis].begin())
238 // Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
239 // else
240 // Log() << Verbose(0) << i << ": " << *runner->second.second;
241 // i++;
242 // }
243 // Log() << Verbose(0) << endl;
244 // }
245 // 3c. throw out points whose distance is less than the mean of left and right neighbours
246 bool flag = false;
247 DoLog(1) && (Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl);
248 do { // do as long as we still throw one out per round
249 flag = false;
250 Boundaries::iterator left = BoundaryPoints[axis].begin();
251 Boundaries::iterator right = BoundaryPoints[axis].begin();
252 Boundaries::iterator runner = BoundaryPoints[axis].begin();
253 bool LoopOnceDone = false;
254 while (!LoopOnceDone) {
255 runner = right;
256 right++;
257 // set neighbours correctly
258 if (runner == BoundaryPoints[axis].begin()) {
259 left = BoundaryPoints[axis].end();
260 } else {
261 left = runner;
262 }
263 left--;
264 if (right == BoundaryPoints[axis].end()) {
265 right = BoundaryPoints[axis].begin();
266 LoopOnceDone = true;
267 }
268 // check distance
269
270 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
271 {
272 Vector SideA, SideB, SideC, SideH;
273 SideA = left->second.second->getPosition() - (*MolCenter);
274 SideA.ProjectOntoPlane(AxisVector);
275 // Log() << Verbose(1) << "SideA: " << SideA << endl;
276
277 SideB = right->second.second->getPosition() -(*MolCenter);
278 SideB.ProjectOntoPlane(AxisVector);
279 // Log() << Verbose(1) << "SideB: " << SideB << endl;
280
281 SideC = left->second.second->getPosition() - right->second.second->getPosition();
282 SideC.ProjectOntoPlane(AxisVector);
283 // Log() << Verbose(1) << "SideC: " << SideC << endl;
284
285 SideH = runner->second.second->getPosition() -(*MolCenter);
286 SideH.ProjectOntoPlane(AxisVector);
287 // Log() << Verbose(1) << "SideH: " << SideH << endl;
288
289 // calculate each length
290 const double a = SideA.Norm();
291 //const double b = SideB.Norm();
292 //const double c = SideC.Norm();
293 const double h = SideH.Norm();
294 // calculate the angles
295 const double alpha = SideA.Angle(SideH);
296 const double beta = SideA.Angle(SideC);
297 const double gamma = SideB.Angle(SideH);
298 const double delta = SideC.Angle(SideH);
299 const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
300 //Log() << Verbose(1) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
301 DoLog(1) && (Log() << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl);
302 if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
303 // throw out point
304 DoLog(1) && (Log() << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl);
305 BoundaryPoints[axis].erase(runner);
306 runner = right;
307 flag = true;
308 }
309 }
310 }
311 } while (flag);
312 }
313 delete(MolCenter);
314 return BoundaryPoints;
315};
316
317/** Tesselates the convex boundary by finding all boundary points.
318 * \param *out output stream for debugging
319 * \param *mol molecule structure with Atom's and Bond's.
320 * \param *BoundaryPts set of boundary points to use or NULL
321 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
322 * \param *LCList atoms in LinkedCell list
323 * \param *filename filename prefix for output of vertex data
324 * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
325 */
326void FindConvexBorder(const molecule* mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
327{
328 Info FunctionInfo(__func__);
329 bool BoundaryFreeFlag = false;
330 Boundaries *BoundaryPoints = NULL;
331
332 if (TesselStruct != NULL) // free if allocated
333 delete(TesselStruct);
334 TesselStruct = new class Tesselation;
335
336 // 1. Find all points on the boundary
337 if (BoundaryPts == NULL) {
338 BoundaryFreeFlag = true;
339 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
340 } else {
341 BoundaryPoints = BoundaryPts;
342 DoLog(0) && (Log() << Verbose(0) << "Using given boundary points set." << endl);
343 }
344
345// printing all inserted for debugging
346 for (int axis=0; axis < NDIM; axis++) {
347 DoLog(1) && (Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl);
348 int i=0;
349 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
350 if (runner != BoundaryPoints[axis].begin())
351 DoLog(0) && (Log() << Verbose(0) << ", " << i << ": " << *runner->second.second);
352 else
353 DoLog(0) && (Log() << Verbose(0) << i << ": " << *runner->second.second);
354 i++;
355 }
356 DoLog(0) && (Log() << Verbose(0) << endl);
357 }
358
359 // 2. fill the boundary point list
360 for (int axis = 0; axis < NDIM; axis++)
361 for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
362 if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
363 DoLog(2) && (Log()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present." << endl);
364
365 DoLog(0) && (Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl);
366 // now we have the whole set of edge points in the BoundaryList
367
368 // listing for debugging
369 // Log() << Verbose(1) << "Listing PointsOnBoundary:";
370 // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
371 // Log() << Verbose(0) << " " << *runner->second;
372 // }
373 // Log() << Verbose(0) << endl;
374
375 // 3a. guess starting triangle
376 TesselStruct->GuessStartingTriangle();
377
378 // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
379 TesselStruct->TesselateOnBoundary(mol);
380
381 // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
382 if (!TesselStruct->InsertStraddlingPoints(mol, LCList))
383 DoeLog(1) && (eLog()<< Verbose(1) << "Insertion of straddling points failed!" << endl);
384
385 DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
386
387 // 4. Store triangles in tecplot file
388 StoreTrianglesinFile(mol, TesselStruct, filename, "_intermed");
389
390 // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
391 bool AllConvex = true;
392 class BoundaryLineSet *line = NULL;
393 do {
394 AllConvex = true;
395 for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
396 line = LineRunner->second;
397 DoLog(1) && (Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl);
398 if (!line->CheckConvexityCriterion()) {
399 DoLog(1) && (Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl);
400
401 // flip the line
402 if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
403 DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
404 else {
405 TesselStruct->FlipBaseline(line);
406 DoLog(1) && (Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl);
407 LineRunner = TesselStruct->LinesOnBoundary.begin(); // LineRunner may have been erase if line was deleted from LinesOnBoundary
408 }
409 }
410 }
411 } while (!AllConvex);
412
413 // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
414// if (!TesselStruct->CorrectConcaveTesselPoints(out))
415// Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;
416
417 DoLog(0) && (Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl);
418
419 // 4. Store triangles in tecplot file
420 StoreTrianglesinFile(mol, TesselStruct, filename, "");
421
422 // free reference lists
423 if (BoundaryFreeFlag)
424 delete[] (BoundaryPoints);
425};
426
427/** For testing removes one boundary point after another to check for leaks.
428 * \param *out output stream for debugging
429 * \param *TesselStruct Tesselation containing envelope with boundary points
430 * \param *mol molecule
431 * \param *filename name of file
432 * \return true - all removed, false - something went wrong
433 */
434bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
435{
436 Info FunctionInfo(__func__);
437 int i=0;
438 char number[MAXSTRINGSIZE];
439
440 if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
441 DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
442 return false;
443 }
444
445 PointMap::iterator PointRunner;
446 while (!TesselStruct->PointsOnBoundary.empty()) {
447 DoLog(1) && (Log() << Verbose(1) << "Remaining points are: ");
448 for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
449 DoLog(0) && (Log() << Verbose(0) << *(PointSprinter->second) << "\t");
450 DoLog(0) && (Log() << Verbose(0) << endl);
451
452 PointRunner = TesselStruct->PointsOnBoundary.begin();
453 // remove point
454 TesselStruct->RemovePointFromTesselatedSurface(PointRunner->second);
455
456 // store envelope
457 sprintf(number, "-%04d", i++);
458 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, number);
459 }
460
461 return true;
462};
463
464/** Creates a convex envelope from a given non-convex one.
465 * -# First step, remove concave spots, i.e. singular "dents"
466 * -# We go through all PointsOnBoundary.
467 * -# We CheckConvexityCriterion() for all its lines.
468 * -# If all its lines are concave, it cannot be on the convex envelope.
469 * -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
470 * -# We calculate the additional volume.
471 * -# We go over all lines until none yields a concavity anymore.
472 * -# Second step, remove concave lines, i.e. line-shape "dents"
473 * -# We go through all LinesOnBoundary
474 * -# We CheckConvexityCriterion()
475 * -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
476 * -# We CheckConvexityCriterion(),
477 * -# if it's concave, we continue
478 * -# if not, we mark an error and stop
479 * Note: This routine - for free - calculates the difference in volume between convex and
480 * non-convex envelope, as the former is easy to calculate - VolumeOfConvexEnvelope() - it
481 * can be used to compute volumes of arbitrary shapes.
482 * \param *out output stream for debugging
483 * \param *TesselStruct non-convex envelope, is changed in return!
484 * \param *mol molecule
485 * \param *filename name of file
486 * \return volume difference between the non- and the created convex envelope
487 */
488double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
489{
490 Info FunctionInfo(__func__);
491 double volume = 0;
492 class BoundaryPointSet *point = NULL;
493 class BoundaryLineSet *line = NULL;
494 bool Concavity = false;
495 char dummy[MAXSTRINGSIZE];
496 PointMap::iterator PointRunner;
497 PointMap::iterator PointAdvance;
498 LineMap::iterator LineRunner;
499 LineMap::iterator LineAdvance;
500 TriangleMap::iterator TriangleRunner;
501 TriangleMap::iterator TriangleAdvance;
502 int run = 0;
503
504 // check whether there is something to work on
505 if (TesselStruct == NULL) {
506 DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
507 return volume;
508 }
509
510 // First step: RemovePointFromTesselatedSurface
511 do {
512 Concavity = false;
513 sprintf(dummy, "-first-%d", run);
514 //CalculateConcavityPerBoundaryPoint(TesselStruct);
515 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
516
517 PointRunner = TesselStruct->PointsOnBoundary.begin();
518 PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
519 while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
520 PointAdvance++;
521 point = PointRunner->second;
522 DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
523 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
524 line = LineRunner->second;
525 DoLog(1) && (Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl);
526 if (!line->CheckConvexityCriterion()) {
527 // remove the point if needed
528 DoLog(1) && (Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl);
529 volume += TesselStruct->RemovePointFromTesselatedSurface(point);
530 sprintf(dummy, "-first-%d", ++run);
531 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
532 Concavity = true;
533 break;
534 }
535 }
536 PointRunner = PointAdvance;
537 }
538
539 sprintf(dummy, "-second-%d", run);
540 //CalculateConcavityPerBoundaryPoint(TesselStruct);
541 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
542
543 // second step: PickFarthestofTwoBaselines
544 LineRunner = TesselStruct->LinesOnBoundary.begin();
545 LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed
546 while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
547 LineAdvance++;
548 line = LineRunner->second;
549 DoLog(1) && (Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl);
550 // take highest of both lines
551 if (TesselStruct->IsConvexRectangle(line) == NULL) {
552 const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
553 volume += tmp;
554 if (tmp != 0.) {
555 TesselStruct->FlipBaseline(line);
556 Concavity = true;
557 }
558 }
559 LineRunner = LineAdvance;
560 }
561 run++;
562 } while (Concavity);
563 //CalculateConcavityPerBoundaryPoint(TesselStruct);
564 //StoreTrianglesinFile(mol, filename, "-third");
565
566 // third step: IsConvexRectangle
567// LineRunner = TesselStruct->LinesOnBoundary.begin();
568// LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed
569// while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
570// LineAdvance++;
571// line = LineRunner->second;
572// Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
573// //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
574// //Log() << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
575// if (!line->CheckConvexityCriterion(out)) {
576// Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
577//
578// // take highest of both lines
579// point = TesselStruct->IsConvexRectangle(line);
580// if (point != NULL)
581// volume += TesselStruct->RemovePointFromTesselatedSurface(point);
582// }
583// LineRunner = LineAdvance;
584// }
585
586 CalculateConcavityPerBoundaryPoint(TesselStruct);
587 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
588
589 // end
590 DoLog(0) && (Log() << Verbose(0) << "Volume is " << volume << "." << endl);
591 return volume;
592};
593
594
595/** Determines the volume of a cluster.
596 * Determines first the convex envelope, then tesselates it and calculates its volume.
597 * \param *out output stream for debugging
598 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
599 * \param *configuration needed for path to store convex envelope file
600 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
601 */
602double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration)
603{
604 Info FunctionInfo(__func__);
605 bool IsAngstroem = configuration->GetIsAngstroem();
606 double volume = 0.;
607 Vector x;
608 Vector y;
609
610 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
611 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
612 { // go through every triangle, calculate volume of its pyramid with CoG as peak
613 x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
614 y = runner->second->getEndpoint(0) - runner->second->getEndpoint(2);
615 const double a = x.Norm();
616 const double b = y.Norm();
617 const double c = runner->second->getEndpoint(2).distance(runner->second->getEndpoint(1));
618 const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
619 x = runner->second->getPlane().getNormal();
620 x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
621 const double h = x.Norm(); // distance of CoG to triangle
622 const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
623 Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
624 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
625 << h << " and the volume is " << PyramidVolume << " "
626 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
627 volume += PyramidVolume;
628 }
629 Log() << Verbose(0) << "RESULT: The summed volume is " << setprecision(6)
630 << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
631 << endl;
632
633 return volume;
634};
635
636/** Stores triangles to file.
637 * \param *out output stream for debugging
638 * \param *mol molecule with atoms and bonds
639 * \param *TesselStruct Tesselation with boundary triangles
640 * \param *filename prefix of filename
641 * \param *extraSuffix intermediate suffix
642 */
643void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
644{
645 Info FunctionInfo(__func__);
646 // 4. Store triangles in tecplot file
647 if (filename != NULL) {
648 if (DoTecplotOutput) {
649 string OutputName(filename);
650 OutputName.append(extraSuffix);
651 OutputName.append(TecplotSuffix);
652 ofstream *tecplot = new ofstream(OutputName.c_str());
653 WriteTecplotFile(tecplot, TesselStruct, mol, -1);
654 tecplot->close();
655 delete(tecplot);
656 }
657 if (DoRaster3DOutput) {
658 string OutputName(filename);
659 OutputName.append(extraSuffix);
660 OutputName.append(Raster3DSuffix);
661 ofstream *rasterplot = new ofstream(OutputName.c_str());
662 WriteRaster3dFile(rasterplot, TesselStruct, mol);
663 rasterplot->close();
664 delete(rasterplot);
665 }
666 }
667};
668
669/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
670 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
671 * TODO: Here, we need a VolumeOfGeneralEnvelope (i.e. non-convex one)
672 * \param *out output stream for debugging
673 * \param *configuration needed for path to store convex envelope file
674 * \param *mol molecule structure representing the cluster
675 * \param *&TesselStruct Tesselation structure with triangles on return
676 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
677 * \param celldensity desired average density in final cell
678 */
679void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
680{
681 Info FunctionInfo(__func__);
682 bool IsAngstroem = true;
683 double *GreatestDiameter = NULL;
684 Boundaries *BoundaryPoints = NULL;
685 class Tesselation *TesselStruct = NULL;
686 Vector BoxLengths;
687 int repetition[NDIM] = { 1, 1, 1 };
688 int TotalNoClusters = 1;
689 double totalmass = 0.;
690 double clustervolume = 0.;
691 double cellvolume = 0.;
692
693 // transform to PAS by Action
694 Vector MainAxis(0.,0.,1.);
695 MoleculeRotateToPrincipalAxisSystem(MainAxis);
696
697 IsAngstroem = configuration->GetIsAngstroem();
698 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
699 GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
700 LinkedCell *LCList = new LinkedCell(*mol, 10.);
701 FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell *&)LCList, NULL);
702 delete (LCList);
703 delete[] BoundaryPoints;
704
705
706 // some preparations beforehand
707 if (ClusterVolume == 0)
708 clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
709 else
710 clustervolume = ClusterVolume;
711
712 delete TesselStruct;
713
714 for (int i = 0; i < NDIM; i++)
715 TotalNoClusters *= repetition[i];
716
717 // sum up the atomic masses
718 for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
719 totalmass += (*iter)->getType()->getMass();
720 }
721 DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
722 DoLog(0) && (Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
723
724 // solve cubic polynomial
725 DoLog(1) && (Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl);
726 if (IsAngstroem)
727 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
728 else
729 cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
730 DoLog(1) && (Log() << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
731
732 double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
733 DoLog(1) && (Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
734 if (minimumvolume > cellvolume) {
735 DoeLog(1) && (eLog()<< Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl);
736 DoLog(0) && (Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl);
737 for (int i = 0; i < NDIM; i++)
738 BoxLengths[i] = GreatestDiameter[i];
739 mol->CenterEdge(&BoxLengths);
740 } else {
741 BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
742 BoxLengths[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
743 BoxLengths[2] = minimumvolume - cellvolume;
744 double x0 = 0.;
745 double x1 = 0.;
746 double x2 = 0.;
747 if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
748 DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl);
749 else {
750 DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl);
751 x0 = x2; // sorted in ascending order
752 }
753
754 cellvolume = 1.;
755 for (int i = 0; i < NDIM; i++) {
756 BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
757 cellvolume *= BoxLengths[i];
758 }
759
760 // set new box dimensions
761 DoLog(0) && (Log() << Verbose(0) << "Translating to box with these boundaries." << endl);
762 mol->SetBoxDimension(&BoxLengths);
763 mol->CenterInBox();
764 }
765 delete GreatestDiameter;
766 // update Box of atoms by boundary
767 mol->SetBoxDimension(&BoxLengths);
768 DoLog(0) && (Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl);
769};
770
771
772/** Fills the empty space around other molecules' surface of the simulation box with a filler.
773 * \param *out output stream for debugging
774 * \param *List list of molecules already present in box
775 * \param *TesselStruct contains tesselated surface
776 * \param *filler molecule which the box is to be filled with
777 * \param configuration contains box dimensions
778 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
779 * \param distance[NDIM] distance between filling molecules in each direction
780 * \param boundary length of boundary zone between molecule and filling mollecules
781 * \param epsilon distance to surface which is not filled
782 * \param RandAtomDisplacement maximum distance for random displacement per atom
783 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
784 * \param DoRandomRotation true - do random rotiations, false - don't
785 * \return *mol pointer to new molecule with filled atoms
786 */
787molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
788{
789 Info FunctionInfo(__func__);
790 molecule *Filling = World::getInstance().createMolecule();
791 Vector CurrentPosition;
792 int N[NDIM];
793 int n[NDIM];
794 const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
795 RealSpaceMatrix Rotations;
796 const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
797 Vector AtomTranslations;
798 Vector FillerTranslations;
799 Vector FillerDistance;
800 Vector Inserter;
801 double FillIt = false;
802 bond *Binder = NULL;
803 double phi[NDIM];
804 map<molecule *, Tesselation *> TesselStruct;
805 map<molecule *, LinkedCell *> LCList;
806
807 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
808 if ((*ListRunner)->getAtomCount() > 0) {
809 DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
810 LCList[(*ListRunner)] = new LinkedCell(*(*ListRunner), 10.); // get linked cell list
811 DoLog(1) && (Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl);
812 TesselStruct[(*ListRunner)] = NULL;
813 FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
814 }
815
816 // Center filler at origin
817 filler->CenterEdge(&Inserter);
818 const int FillerCount = filler->getAtomCount();
819 DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
820 for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
821 for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
822 if ((*BondRunner)->leftatom == *AtomRunner)
823 DoLog(2) && (Log() << Verbose(2) << " " << *(*BondRunner) << endl);
824
825 atom * CopyAtoms[FillerCount];
826
827 // calculate filler grid in [0,1]^3
828 FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
829 for(int i=0;i<NDIM;i++)
830 N[i] = (int) ceil(1./FillerDistance[i]);
831 DoLog(1) && (Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl);
832
833 // initialize seed of random number generator to current time
834 srand ( time(NULL) );
835
836 // go over [0,1]^3 filler grid
837 for (n[0] = 0; n[0] < N[0]; n[0]++)
838 for (n[1] = 0; n[1] < N[1]; n[1]++)
839 for (n[2] = 0; n[2] < N[2]; n[2]++) {
840 // calculate position of current grid vector in untransformed box
841 CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
842 // create molecule random translation vector ...
843 for (int i=0;i<NDIM;i++)
844 FillerTranslations[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
845 DoLog(2) && (Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl);
846
847 // go through all atoms
848 for (int i=0;i<FillerCount;i++)
849 CopyAtoms[i] = NULL;
850
851 // have same rotation angles for all molecule's atoms
852 if (DoRandomRotation)
853 for (int i=0;i<NDIM;i++)
854 phi[i] = (rand()/RAND_MAX)*(2.*M_PI);
855
856 for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){
857
858 // create atomic random translation vector ...
859 for (int i=0;i<NDIM;i++)
860 AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
861
862 // ... and rotation matrix
863 if (DoRandomRotation) {
864 Rotations.set(0,0, cos(phi[0]) *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
865 Rotations.set(0,1, sin(phi[0]) *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
866 Rotations.set(0,2, cos(phi[1])*sin(phi[2]) );
867 Rotations.set(1,0, -sin(phi[0])*cos(phi[1]) );
868 Rotations.set(1,1, cos(phi[0])*cos(phi[1]) );
869 Rotations.set(1,2, sin(phi[1]) );
870 Rotations.set(2,0, -cos(phi[0]) *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
871 Rotations.set(2,1, -sin(phi[0]) *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
872 Rotations.set(2,2, cos(phi[1])*cos(phi[2]) );
873 }
874
875 // ... and put at new position
876 Inserter = (*iter)->getPosition();
877 if (DoRandomRotation)
878 Inserter *= Rotations;
879 Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
880
881 // check whether inserter is inside box
882 Inserter *= MInverse;
883 FillIt = true;
884 for (int i=0;i<NDIM;i++)
885 FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
886 Inserter *= M;
887
888 // Check whether point is in- or outside
889 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
890 // get linked cell list
891 if (TesselStruct[(*ListRunner)] != NULL) {
892 const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
893 FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
894 }
895 }
896 // insert into Filling
897 if (FillIt) {
898 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl);
899 // copy atom ...
900 CopyAtoms[(*iter)->nr] = (*iter)->clone();
901 (*CopyAtoms[(*iter)->nr]).setPosition(Inserter);
902 Filling->AddAtom(CopyAtoms[(*iter)->nr]);
903 DoLog(1) && (Log() << Verbose(1) << "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->getPosition()) << "." << endl);
904 } else {
905 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);
906 CopyAtoms[(*iter)->nr] = NULL;
907 continue;
908 }
909 }
910 // go through all bonds and add as well
911 for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
912 for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
913 if ((*BondRunner)->leftatom == *AtomRunner) {
914 Binder = (*BondRunner);
915 if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
916 Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
917 Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
918 }
919 }
920 }
921 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
922 delete LCList[*ListRunner];
923 delete TesselStruct[(*ListRunner)];
924 }
925
926 return Filling;
927};
928
929/** Rotates given molecule \a Filling and moves its atoms according to given
930 * \a RandomAtomDisplacement.
931 *
932 * Note that for rotation to be sensible, the molecule should be centered at
933 * the origin. This is not done here!
934 *
935 * \param &Filling molecule whose atoms to displace
936 * \param RandomAtomDisplacement magnitude of random displacement
937 * \param &Rotations 3D rotation matrix (or unity if no rotation desired)
938 */
939void RandomizeMoleculePositions(
940 molecule *&Filling,
941 double RandomAtomDisplacement,
942 RealSpaceMatrix &Rotations
943 )
944{
945 Vector AtomTranslations;
946 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
947 Vector temp = (*miter)->getPosition();
948 temp *= Rotations;
949 (*miter)->setPosition(temp);
950 // create atomic random translation vector ...
951 for (int i=0;i<NDIM;i++)
952 AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
953 (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
954 }
955}
956
957/** Removes all atoms of a molecule outside.
958 *
959 * If the molecule is empty, it is removed as well.
960 *
961 * @param Filling molecule whose atoms to check, removed if eventually left
962 * empty.
963 * @return true - atoms had to be removed, false - no atoms have been removed
964 */
965bool RemoveAtomsOutsideDomain(molecule *&Filling)
966{
967 bool status = false;
968 Box &Domain = World::getInstance().getDomain();
969 // check if all is still inside domain
970 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
971 // check whether each atom is inside box
972 if (!Domain.isInside((*miter)->getPosition())) {
973 status = true;
974 atom *Walker = *miter;
975 ++miter;
976 World::getInstance().destroyAtom(Walker);
977 } else {
978 ++miter;
979 }
980 }
981 if (Filling->empty()) {
982 DoLog(0) && (Log() << Verbose(0) << "Removing molecule " << Filling->getName() << ", all atoms have been removed." << std::endl);
983 World::getInstance().destroyMolecule(Filling);
984 }
985 return status;
986}
987
988/** Checks whether there are no atoms inside a sphere around \a CurrentPosition
989 * except those atoms present in \a *filler.
990 * If filler is NULL, then we just call LinkedCell::GetPointsInsideSphere() and
991 * check whether the return list is empty.
992 * @param *filler
993 * @param boundary
994 * @param CurrentPosition
995 */
996bool isSpaceAroundPointVoid(
997 LinkedCell *LC,
998 molecule *filler,
999 const double boundary,
1000 Vector &CurrentPosition)
1001{
1002 size_t compareTo = 0;
1003 LinkedCell::LinkedNodes* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition);
1004 if (filler != NULL) {
1005 for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin();
1006 iter != liste->end();
1007 ++iter) {
1008 for (molecule::iterator miter = filler->begin();
1009 miter != filler->end();
1010 ++miter) {
1011 if (*iter == *miter)
1012 ++compareTo;
1013 }
1014 }
1015 }
1016 const bool result = (liste->size() == compareTo);
1017 if (!result) {
1018 DoLog(0) && (Log() << Verbose(0) << "Skipping because of the following atoms:" << std::endl);
1019 for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin();
1020 iter != liste->end();
1021 ++iter) {
1022 DoLog(0) && (Log() << Verbose(0) << **iter << std::endl);
1023 }
1024 }
1025 delete(liste);
1026 return result;
1027}
1028
1029/** Fills the empty space of the simulation box with water.
1030 * \param *filler molecule which the box is to be filled with
1031 * \param configuration contains box dimensions
1032 * \param distance[NDIM] distance between filling molecules in each direction
1033 * \param boundary length of boundary zone between molecule and filling molecules
1034 * \param RandAtomDisplacement maximum distance for random displacement per atom
1035 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
1036 * \param MinDistance minimum distance to boundary of domain and present molecules
1037 * \param DoRandomRotation true - do random rotations, false - don't
1038 */
1039void FillVoidWithMolecule(
1040 molecule *&filler,
1041 config &configuration,
1042 const double distance[NDIM],
1043 const double boundary,
1044 const double RandomAtomDisplacement,
1045 const double RandomMolDisplacement,
1046 const double MinDistance,
1047 const bool DoRandomRotation
1048 )
1049{
1050 Info FunctionInfo(__func__);
1051 molecule *Filling = NULL;
1052 Vector CurrentPosition;
1053 int N[NDIM];
1054 int n[NDIM];
1055 const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
1056 RealSpaceMatrix Rotations;
1057 const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
1058 Vector FillerTranslations;
1059 Vector FillerDistance;
1060 Vector Inserter;
1061 double FillIt = false;
1062 Vector firstInserter;
1063 bool firstInsertion = true;
1064 const Box &Domain = World::getInstance().getDomain();
1065 map<molecule *, LinkedCell *> LCList;
1066 std::vector<molecule *> List = World::getInstance().getAllMolecules();
1067 MoleculeListClass *MolList = World::getInstance().getMolecules();
1068
1069 for (std::vector<molecule *>::iterator ListRunner = List.begin(); ListRunner != List.end(); ListRunner++)
1070 if ((*ListRunner)->getAtomCount() > 0) {
1071 DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
1072 LCList[(*ListRunner)] = new LinkedCell(*(*ListRunner), 10.); // get linked cell list
1073 }
1074
1075 // Center filler at its center of gravity
1076 Vector *gravity = filler->DetermineCenterOfGravity();
1077 filler->CenterAtVector(gravity);
1078 delete gravity;
1079 //const int FillerCount = filler->getAtomCount();
1080 DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
1081 for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
1082 for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
1083 if ((*BondRunner)->leftatom == *AtomRunner)
1084 DoLog(2) && (Log() << Verbose(2) << " " << *(*BondRunner) << endl);
1085
1086 // calculate filler grid in [0,1]^3
1087 FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
1088 for(int i=0;i<NDIM;i++)
1089 N[i] = (int) ceil(1./FillerDistance[i]);
1090 DoLog(1) && (Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl);
1091
1092 // initialize seed of random number generator to current time
1093 srand ( time(NULL) );
1094
1095 // go over [0,1]^3 filler grid
1096 for (n[0] = 0; n[0] < N[0]; n[0]++)
1097 for (n[1] = 0; n[1] < N[1]; n[1]++)
1098 for (n[2] = 0; n[2] < N[2]; n[2]++) {
1099 // calculate position of current grid vector in untransformed box
1100 CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
1101 // create molecule random translation vector ...
1102 for (int i=0;i<NDIM;i++)
1103 FillerTranslations[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
1104 DoLog(2) && (Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl);
1105
1106 // ... and rotation matrix
1107 if (DoRandomRotation)
1108 Rotations.setRandomRotation();
1109 else
1110 Rotations.setIdentity();
1111
1112
1113 // Check whether there is anything too close by and whether atom is outside of domain
1114 FillIt = true;
1115 for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
1116 FillIt = FillIt && isSpaceAroundPointVoid(
1117 ListRunner->second,
1118 (firstInsertion ? filler : NULL),
1119 boundary,
1120 CurrentPosition);
1121 FillIt = FillIt && (Domain.isInside(CurrentPosition))
1122 && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
1123 if (!FillIt)
1124 break;
1125 }
1126
1127 // insert into Filling
1128 if (FillIt) {
1129 Inserter = CurrentPosition + FillerTranslations;
1130 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is void point." << endl);
1131 // fill!
1132 Filling = filler->CopyMolecule();
1133 RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations);
1134 // translation
1135 Filling->Translate(&Inserter);
1136 // remove out-of-bounds atoms
1137 const bool status = RemoveAtomsOutsideDomain(Filling);
1138 if ((firstInsertion) && (!status)) { // no atom has been removed
1139 // remove copied atoms and molecule again
1140 Filling->removeAtomsinMolecule();
1141 World::getInstance().destroyMolecule(Filling);
1142 // and mark is final filler position
1143 Filling = filler;
1144 firstInsertion = false;
1145 firstInserter = Inserter;
1146 } else {
1147 // TODO: Remove when World has no MoleculeListClass anymore
1148 MolList->insert(Filling);
1149 }
1150 } else {
1151 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance." << endl);
1152 continue;
1153 }
1154 }
1155
1156 // have we inserted any molecules?
1157 if (firstInsertion) {
1158 // If not remove filler
1159 for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
1160 atom *Walker = *miter;
1161 filler->erase(Walker);
1162 World::getInstance().destroyAtom(Walker);
1163 }
1164 World::getInstance().destroyMolecule(filler);
1165 } else {
1166 // otherwise translate and randomize to final position
1167 if (DoRandomRotation)
1168 Rotations.setRandomRotation();
1169 else
1170 Rotations.setIdentity();
1171 RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations);
1172 // translation
1173 filler->Translate(&firstInserter);
1174 // remove out-of-bounds atoms
1175 RemoveAtomsOutsideDomain(filler);
1176 }
1177
1178 DoLog(0) && (Log() << Verbose(0) << MolList->ListOfMolecules.size() << " molecules have been inserted." << std::endl);
1179
1180 for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
1181 delete ListRunner->second;
1182 LCList.erase(ListRunner);
1183 }
1184};
1185
1186/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
1187 * \param *out output stream for debugging
1188 * \param *mol molecule structure with Atom's and Bond's
1189 * \param *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
1190 * \param *&LCList atoms in LinkedCell list
1191 * \param RADIUS radius of the virtual sphere
1192 * \param *filename filename prefix for output of vertex data
1193 * \return true - tesselation successful, false - tesselation failed
1194 */
1195bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
1196{
1197 Info FunctionInfo(__func__);
1198 bool freeLC = false;
1199 bool status = false;
1200 CandidateForTesselation *baseline = NULL;
1201 bool OneLoopWithoutSuccessFlag = true; // marks whether we went once through all baselines without finding any without two triangles
1202 bool TesselationFailFlag = false;
1203
1204 mol->getAtomCount();
1205
1206 if (TesselStruct == NULL) {
1207 DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
1208 TesselStruct= new Tesselation;
1209 } else {
1210 delete(TesselStruct);
1211 DoLog(1) && (Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl);
1212 TesselStruct = new Tesselation;
1213 }
1214
1215 // initialise Linked Cell
1216 if (LCList == NULL) {
1217 LCList = new LinkedCell(*mol, 2.*RADIUS);
1218 freeLC = true;
1219 }
1220
1221 // 1. get starting triangle
1222 if (!TesselStruct->FindStartingTriangle(RADIUS, LCList)) {
1223 DoeLog(0) && (eLog() << Verbose(0) << "No valid starting triangle found." << endl);
1224 //performCriticalExit();
1225 }
1226 if (filename != NULL) {
1227 if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
1228 TesselStruct->Output(filename, mol);
1229 }
1230 }
1231
1232 // 2. expand from there
1233 while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
1234 (cerr << "There are " << TesselStruct->TrianglesOnBoundary.size() << " triangles and " << TesselStruct->OpenLines.size() << " open lines to scan for candidates." << endl);
1235 // 2a. print OpenLines without candidates
1236 DoLog(1) && (Log() << Verbose(1) << "There are the following open lines to scan for a candidates:" << endl);
1237 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
1238 if (Runner->second->pointlist.empty())
1239 DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
1240
1241 // 2b. find best candidate for each OpenLine
1242 TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);
1243
1244 // 2c. print OpenLines with candidates again
1245 DoLog(1) && (Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl);
1246 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
1247 DoLog(1) && (Log() << Verbose(1) << " " << *(Runner->second) << endl);
1248
1249 // 2d. search for smallest ShortestAngle among all candidates
1250 double ShortestAngle = 4.*M_PI;
1251 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
1252 if (Runner->second->ShortestAngle < ShortestAngle) {
1253 baseline = Runner->second;
1254 ShortestAngle = baseline->ShortestAngle;
1255 DoLog(1) && (Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl);
1256 }
1257 }
1258 // 2e. if we found one, add candidate
1259 if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
1260 OneLoopWithoutSuccessFlag = false;
1261 else {
1262 TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
1263 }
1264
1265 // 2f. write temporary envelope
1266 if (filename != NULL) {
1267 if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
1268 TesselStruct->Output(filename, mol);
1269 }
1270 }
1271 }
1272// // check envelope for consistency
1273// status = CheckListOfBaselines(TesselStruct);
1274//
1275// // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
1276// //->InsertStraddlingPoints(mol, LCList);
1277// for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
1278// class TesselPoint *Runner = NULL;
1279// Runner = *iter;
1280// Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
1281// if (!->IsInnerPoint(Runner, LCList)) {
1282// Log() << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
1283// ->AddBoundaryPointByDegeneratedTriangle(Runner, LCList);
1284// } else {
1285// Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
1286// }
1287// }
1288
1289// // Purges surplus triangles.
1290// TesselStruct->RemoveDegeneratedTriangles();
1291//
1292// // check envelope for consistency
1293// status = CheckListOfBaselines(TesselStruct);
1294
1295 cout << "before correction" << endl;
1296
1297 // store before correction
1298 StoreTrianglesinFile(mol, TesselStruct, filename, "");
1299
1300// // correct degenerated polygons
1301// TesselStruct->CorrectAllDegeneratedPolygons();
1302//
1303// // check envelope for consistency
1304// status = CheckListOfBaselines(TesselStruct);
1305
1306 // write final envelope
1307 CalculateConcavityPerBoundaryPoint(TesselStruct);
1308 cout << "after correction" << endl;
1309 StoreTrianglesinFile(mol, TesselStruct, filename, "");
1310
1311 if (freeLC)
1312 delete(LCList);
1313
1314 return status;
1315};
1316
1317
1318/** Finds a hole of sufficient size in \a *mols to embed \a *srcmol into it.
1319 * \param *out output stream for debugging
1320 * \param *mols molecules in the domain to embed in between
1321 * \param *srcmol embedding molecule
1322 * \return *Vector new center of \a *srcmol for embedding relative to \a this
1323 */
1324Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
1325{
1326 Info FunctionInfo(__func__);
1327 Vector *Center = new Vector;
1328 Center->Zero();
1329 // calculate volume/shape of \a *srcmol
1330
1331 // find embedding holes
1332
1333 // if more than one, let user choose
1334
1335 // return embedding center
1336 return Center;
1337};
1338
Note: See TracBrowser for help on using the repository browser.