source: src/boundary.cpp@ b8c9d8

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

Don't exit when no starting triangle has been found.

  • The issue is that this is not a bug, but just a limitation of the so-far employed algorithm for obtaining a starting triangle.
  • If we exit, no .dat file is written, not even an empty one, indicating a perpetual running of the routine, which it does not.

Signed-off-by: Frederik Heber <heber@…>

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