| 1 | /** \file boundary.hpp
|
|---|
| 2 | *
|
|---|
| 3 | * Implementations and super-function for envelopes
|
|---|
| 4 | */
|
|---|
| 5 |
|
|---|
| 6 |
|
|---|
| 7 | #include "boundary.hpp"
|
|---|
| 8 |
|
|---|
| 9 | #include<gsl/gsl_poly.h>
|
|---|
| 10 |
|
|---|
| 11 | // ========================================== F U N C T I O N S =================================
|
|---|
| 12 |
|
|---|
| 13 |
|
|---|
| 14 | /** Determines greatest diameters of a cluster defined by its convex envelope.
|
|---|
| 15 | * Looks at lines parallel to one axis and where they intersect on the projected planes
|
|---|
| 16 | * \param *out output stream for debugging
|
|---|
| 17 | * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
|
|---|
| 18 | * \param *mol molecule structure representing the cluster
|
|---|
| 19 | * \param IsAngstroem whether we have angstroem or atomic units
|
|---|
| 20 | * \return NDIM array of the diameters
|
|---|
| 21 | */
|
|---|
| 22 | double *
|
|---|
| 23 | GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
|
|---|
| 24 | bool IsAngstroem)
|
|---|
| 25 | {
|
|---|
| 26 | // get points on boundary of NULL was given as parameter
|
|---|
| 27 | bool BoundaryFreeFlag = false;
|
|---|
| 28 | Boundaries *BoundaryPoints = BoundaryPtr;
|
|---|
| 29 | if (BoundaryPoints == NULL) {
|
|---|
| 30 | BoundaryFreeFlag = true;
|
|---|
| 31 | BoundaryPoints = GetBoundaryPoints(out, mol);
|
|---|
| 32 | } else {
|
|---|
| 33 | *out << Verbose(1) << "Using given boundary points set." << endl;
|
|---|
| 34 | }
|
|---|
| 35 | // determine biggest "diameter" of cluster for each axis
|
|---|
| 36 | Boundaries::iterator Neighbour, OtherNeighbour;
|
|---|
| 37 | double *GreatestDiameter = new double[NDIM];
|
|---|
| 38 | for (int i = 0; i < NDIM; i++)
|
|---|
| 39 | GreatestDiameter[i] = 0.;
|
|---|
| 40 | double OldComponent, tmp, w1, w2;
|
|---|
| 41 | Vector DistanceVector, OtherVector;
|
|---|
| 42 | int component, Othercomponent;
|
|---|
| 43 | for (int axis = 0; axis < NDIM; axis++)
|
|---|
| 44 | { // regard each projected plane
|
|---|
| 45 | //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
|
|---|
| 46 | for (int j = 0; j < 2; j++)
|
|---|
| 47 | { // and for both axis on the current plane
|
|---|
| 48 | component = (axis + j + 1) % NDIM;
|
|---|
| 49 | Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
|
|---|
| 50 | //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
|
|---|
| 51 | for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
|
|---|
| 52 | != BoundaryPoints[axis].end(); runner++)
|
|---|
| 53 | {
|
|---|
| 54 | //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
|
|---|
| 55 | // seek for the neighbours pair where the Othercomponent sign flips
|
|---|
| 56 | Neighbour = runner;
|
|---|
| 57 | Neighbour++;
|
|---|
| 58 | if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
|
|---|
| 59 | Neighbour = BoundaryPoints[axis].begin();
|
|---|
| 60 | DistanceVector.CopyVector(&runner->second.second->x);
|
|---|
| 61 | DistanceVector.SubtractVector(&Neighbour->second.second->x);
|
|---|
| 62 | do
|
|---|
| 63 | { // seek for neighbour pair where it flips
|
|---|
| 64 | OldComponent = DistanceVector.x[Othercomponent];
|
|---|
| 65 | Neighbour++;
|
|---|
| 66 | if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
|
|---|
| 67 | Neighbour = BoundaryPoints[axis].begin();
|
|---|
| 68 | DistanceVector.CopyVector(&runner->second.second->x);
|
|---|
| 69 | DistanceVector.SubtractVector(&Neighbour->second.second->x);
|
|---|
| 70 | //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
|
|---|
| 71 | }
|
|---|
| 72 | while ((runner != Neighbour) && (fabs(OldComponent / fabs(
|
|---|
| 73 | OldComponent) - DistanceVector.x[Othercomponent] / fabs(
|
|---|
| 74 | DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
|
|---|
| 75 | if (runner != Neighbour)
|
|---|
| 76 | {
|
|---|
| 77 | OtherNeighbour = Neighbour;
|
|---|
| 78 | if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
|
|---|
| 79 | OtherNeighbour = BoundaryPoints[axis].end();
|
|---|
| 80 | OtherNeighbour--;
|
|---|
| 81 | //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
|
|---|
| 82 | // now we have found the pair: Neighbour and OtherNeighbour
|
|---|
| 83 | OtherVector.CopyVector(&runner->second.second->x);
|
|---|
| 84 | OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
|
|---|
| 85 | //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
|
|---|
| 86 | //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
|
|---|
| 87 | // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
|
|---|
| 88 | w1 = fabs(OtherVector.x[Othercomponent]);
|
|---|
| 89 | w2 = fabs(DistanceVector.x[Othercomponent]);
|
|---|
| 90 | tmp = fabs((w1 * DistanceVector.x[component] + w2
|
|---|
| 91 | * OtherVector.x[component]) / (w1 + w2));
|
|---|
| 92 | // mark if it has greater diameter
|
|---|
| 93 | //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
|
|---|
| 94 | GreatestDiameter[component] = (GreatestDiameter[component]
|
|---|
| 95 | > tmp) ? GreatestDiameter[component] : tmp;
|
|---|
| 96 | } //else
|
|---|
| 97 | //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
|
|---|
| 98 | }
|
|---|
| 99 | }
|
|---|
| 100 | }
|
|---|
| 101 | *out << Verbose(0) << "RESULT: The biggest diameters are "
|
|---|
| 102 | << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
|
|---|
| 103 | << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
|
|---|
| 104 | : "atomiclength") << "." << endl;
|
|---|
| 105 |
|
|---|
| 106 | // free reference lists
|
|---|
| 107 | if (BoundaryFreeFlag)
|
|---|
| 108 | delete[] (BoundaryPoints);
|
|---|
| 109 |
|
|---|
| 110 | return GreatestDiameter;
|
|---|
| 111 | }
|
|---|
| 112 | ;
|
|---|
| 113 |
|
|---|
| 114 | /** Creates the objects in a VRML file.
|
|---|
| 115 | * \param *out output stream for debugging
|
|---|
| 116 | * \param *vrmlfile output stream for tecplot data
|
|---|
| 117 | * \param *Tess Tesselation structure with constructed triangles
|
|---|
| 118 | * \param *mol molecule structure with atom positions
|
|---|
| 119 | */
|
|---|
| 120 | void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
|
|---|
| 121 | {
|
|---|
| 122 | atom *Walker = mol->start;
|
|---|
| 123 | bond *Binder = mol->first;
|
|---|
| 124 | int i;
|
|---|
| 125 | Vector *center = mol->DetermineCenterOfAll(out);
|
|---|
| 126 | if (vrmlfile != NULL) {
|
|---|
| 127 | //cout << Verbose(1) << "Writing Raster3D file ... ";
|
|---|
| 128 | *vrmlfile << "#VRML V2.0 utf8" << endl;
|
|---|
| 129 | *vrmlfile << "#Created by molecuilder" << endl;
|
|---|
| 130 | *vrmlfile << "#All atoms as spheres" << endl;
|
|---|
| 131 | while (Walker->next != mol->end) {
|
|---|
| 132 | Walker = Walker->next;
|
|---|
| 133 | *vrmlfile << "Sphere {" << endl << " "; // 2 is sphere type
|
|---|
| 134 | for (i=0;i<NDIM;i++)
|
|---|
| 135 | *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
|
|---|
| 136 | *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
|
|---|
| 137 | }
|
|---|
| 138 |
|
|---|
| 139 | *vrmlfile << "# All bonds as vertices" << endl;
|
|---|
| 140 | while (Binder->next != mol->last) {
|
|---|
| 141 | Binder = Binder->next;
|
|---|
| 142 | *vrmlfile << "3" << endl << " "; // 2 is round-ended cylinder type
|
|---|
| 143 | for (i=0;i<NDIM;i++)
|
|---|
| 144 | *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
|
|---|
| 145 | *vrmlfile << "\t0.03\t";
|
|---|
| 146 | for (i=0;i<NDIM;i++)
|
|---|
| 147 | *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
|
|---|
| 148 | *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
|
|---|
| 149 | }
|
|---|
| 150 |
|
|---|
| 151 | *vrmlfile << "# All tesselation triangles" << endl;
|
|---|
| 152 | for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
|
|---|
| 153 | *vrmlfile << "1" << endl << " "; // 1 is triangle type
|
|---|
| 154 | for (i=0;i<3;i++) { // print each node
|
|---|
| 155 | for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates
|
|---|
| 156 | *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
|
|---|
| 157 | *vrmlfile << "\t";
|
|---|
| 158 | }
|
|---|
| 159 | *vrmlfile << "1. 0. 0." << endl; // red as colour
|
|---|
| 160 | *vrmlfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
|---|
| 161 | }
|
|---|
| 162 | } else {
|
|---|
| 163 | cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
|
|---|
| 164 | }
|
|---|
| 165 | delete(center);
|
|---|
| 166 | };
|
|---|
| 167 |
|
|---|
| 168 | /** Creates the objects in a raster3d file (renderable with a header.r3d).
|
|---|
| 169 | * \param *out output stream for debugging
|
|---|
| 170 | * \param *rasterfile output stream for tecplot data
|
|---|
| 171 | * \param *Tess Tesselation structure with constructed triangles
|
|---|
| 172 | * \param *mol molecule structure with atom positions
|
|---|
| 173 | */
|
|---|
| 174 | void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
|
|---|
| 175 | {
|
|---|
| 176 | atom *Walker = mol->start;
|
|---|
| 177 | bond *Binder = mol->first;
|
|---|
| 178 | int i;
|
|---|
| 179 | Vector *center = mol->DetermineCenterOfAll(out);
|
|---|
| 180 | if (rasterfile != NULL) {
|
|---|
| 181 | //cout << Verbose(1) << "Writing Raster3D file ... ";
|
|---|
| 182 | *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
|
|---|
| 183 | *rasterfile << "@header.r3d" << endl;
|
|---|
| 184 | *rasterfile << "# All atoms as spheres" << endl;
|
|---|
| 185 | while (Walker->next != mol->end) {
|
|---|
| 186 | Walker = Walker->next;
|
|---|
| 187 | *rasterfile << "2" << endl << " "; // 2 is sphere type
|
|---|
| 188 | for (i=0;i<NDIM;i++)
|
|---|
| 189 | *rasterfile << Walker->x.x[i]-center->x[i] << " ";
|
|---|
| 190 | *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
|
|---|
| 191 | }
|
|---|
| 192 |
|
|---|
| 193 | *rasterfile << "# All bonds as vertices" << endl;
|
|---|
| 194 | while (Binder->next != mol->last) {
|
|---|
| 195 | Binder = Binder->next;
|
|---|
| 196 | *rasterfile << "3" << endl << " "; // 2 is round-ended cylinder type
|
|---|
| 197 | for (i=0;i<NDIM;i++)
|
|---|
| 198 | *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
|
|---|
| 199 | *rasterfile << "\t0.03\t";
|
|---|
| 200 | for (i=0;i<NDIM;i++)
|
|---|
| 201 | *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
|
|---|
| 202 | *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
|
|---|
| 203 | }
|
|---|
| 204 |
|
|---|
| 205 | *rasterfile << "# All tesselation triangles" << endl;
|
|---|
| 206 | *rasterfile << "8\n 25. -1. 1. 1. 1. 0.0 0 0 0 2\n SOLID 1.0 0.0 0.0\n BACKFACE 0.3 0.3 1.0 0 0\n";
|
|---|
| 207 | for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
|
|---|
| 208 | *rasterfile << "1" << endl << " "; // 1 is triangle type
|
|---|
| 209 | for (i=0;i<3;i++) { // print each node
|
|---|
| 210 | for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates
|
|---|
| 211 | *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
|
|---|
| 212 | *rasterfile << "\t";
|
|---|
| 213 | }
|
|---|
| 214 | *rasterfile << "1. 0. 0." << endl; // red as colour
|
|---|
| 215 | //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
|---|
| 216 | }
|
|---|
| 217 | *rasterfile << "9\n terminating special property\n";
|
|---|
| 218 | } else {
|
|---|
| 219 | cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
|
|---|
| 220 | }
|
|---|
| 221 | delete(center);
|
|---|
| 222 | };
|
|---|
| 223 |
|
|---|
| 224 | /** This function creates the tecplot file, displaying the tesselation of the hull.
|
|---|
| 225 | * \param *out output stream for debugging
|
|---|
| 226 | * \param *tecplot output stream for tecplot data
|
|---|
| 227 | * \param N arbitrary number to differentiate various zones in the tecplot format
|
|---|
| 228 | */
|
|---|
| 229 | void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
|
|---|
| 230 | {
|
|---|
| 231 | if (tecplot != NULL) {
|
|---|
| 232 | *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
|
|---|
| 233 | *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
|
|---|
| 234 | *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
|
|---|
| 235 | int *LookupList = new int[mol->AtomCount];
|
|---|
| 236 | for (int i = 0; i < mol->AtomCount; i++)
|
|---|
| 237 | LookupList[i] = -1;
|
|---|
| 238 |
|
|---|
| 239 | // print atom coordinates
|
|---|
| 240 | *out << Verbose(2) << "The following triangles were created:";
|
|---|
| 241 | int Counter = 1;
|
|---|
| 242 | TesselPoint *Walker = NULL;
|
|---|
| 243 | for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
|
|---|
| 244 | Walker = target->second->node;
|
|---|
| 245 | LookupList[Walker->nr] = Counter++;
|
|---|
| 246 | *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
|
|---|
| 247 | }
|
|---|
| 248 | *tecplot << endl;
|
|---|
| 249 | // print connectivity
|
|---|
| 250 | for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
|
|---|
| 251 | *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
|
|---|
| 252 | *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
|
|---|
| 253 | }
|
|---|
| 254 | delete[] (LookupList);
|
|---|
| 255 | *out << endl;
|
|---|
| 256 | }
|
|---|
| 257 | }
|
|---|
| 258 |
|
|---|
| 259 |
|
|---|
| 260 | /** Determines the boundary points of a cluster.
|
|---|
| 261 | * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
|
|---|
| 262 | * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
|
|---|
| 263 | * center and first and last point in the triple, it is thrown out.
|
|---|
| 264 | * \param *out output stream for debugging
|
|---|
| 265 | * \param *mol molecule structure representing the cluster
|
|---|
| 266 | */
|
|---|
| 267 | Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
|
|---|
| 268 | {
|
|---|
| 269 | atom *Walker = NULL;
|
|---|
| 270 | PointMap PointsOnBoundary;
|
|---|
| 271 | LineMap LinesOnBoundary;
|
|---|
| 272 | TriangleMap TrianglesOnBoundary;
|
|---|
| 273 | Vector *MolCenter = mol->DetermineCenterOfAll(out);
|
|---|
| 274 | Vector helper;
|
|---|
| 275 |
|
|---|
| 276 | *out << Verbose(1) << "Finding all boundary points." << endl;
|
|---|
| 277 | Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
|
|---|
| 278 | BoundariesTestPair BoundaryTestPair;
|
|---|
| 279 | Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
|
|---|
| 280 | double radius, angle;
|
|---|
| 281 | // 3a. Go through every axis
|
|---|
| 282 | for (int axis = 0; axis < NDIM; axis++) {
|
|---|
| 283 | AxisVector.Zero();
|
|---|
| 284 | AngleReferenceVector.Zero();
|
|---|
| 285 | AngleReferenceNormalVector.Zero();
|
|---|
| 286 | AxisVector.x[axis] = 1.;
|
|---|
| 287 | AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
|
|---|
| 288 | AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
|
|---|
| 289 |
|
|---|
| 290 | *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
|
|---|
| 291 |
|
|---|
| 292 | // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
|
|---|
| 293 | Walker = mol->start;
|
|---|
| 294 | while (Walker->next != mol->end) {
|
|---|
| 295 | Walker = Walker->next;
|
|---|
| 296 | Vector ProjectedVector;
|
|---|
| 297 | ProjectedVector.CopyVector(&Walker->x);
|
|---|
| 298 | ProjectedVector.SubtractVector(MolCenter);
|
|---|
| 299 | ProjectedVector.ProjectOntoPlane(&AxisVector);
|
|---|
| 300 |
|
|---|
| 301 | // correct for negative side
|
|---|
| 302 | radius = ProjectedVector.NormSquared();
|
|---|
| 303 | if (fabs(radius) > MYEPSILON)
|
|---|
| 304 | angle = ProjectedVector.Angle(&AngleReferenceVector);
|
|---|
| 305 | else
|
|---|
| 306 | angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
|
|---|
| 307 |
|
|---|
| 308 | //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
|
|---|
| 309 | if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
|
|---|
| 310 | angle = 2. * M_PI - angle;
|
|---|
| 311 | }
|
|---|
| 312 | *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
|
|---|
| 313 | BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
|
|---|
| 314 | if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
|
|---|
| 315 | *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
|
|---|
| 316 | *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
|
|---|
| 317 | *out << Verbose(2) << "New vector: " << *Walker << endl;
|
|---|
| 318 | double tmp = ProjectedVector.NormSquared();
|
|---|
| 319 | if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
|
|---|
| 320 | BoundaryTestPair.first->second.first = tmp;
|
|---|
| 321 | BoundaryTestPair.first->second.second = Walker;
|
|---|
| 322 | *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
|
|---|
| 323 | } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
|
|---|
| 324 | helper.CopyVector(&Walker->x);
|
|---|
| 325 | helper.SubtractVector(MolCenter);
|
|---|
| 326 | tmp = helper.NormSquared();
|
|---|
| 327 | helper.CopyVector(&BoundaryTestPair.first->second.second->x);
|
|---|
| 328 | helper.SubtractVector(MolCenter);
|
|---|
| 329 | if (helper.NormSquared() < tmp) {
|
|---|
| 330 | BoundaryTestPair.first->second.second = Walker;
|
|---|
| 331 | *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
|
|---|
| 332 | } else {
|
|---|
| 333 | *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
|
|---|
| 334 | }
|
|---|
| 335 | } else {
|
|---|
| 336 | *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
|
|---|
| 337 | }
|
|---|
| 338 | }
|
|---|
| 339 | }
|
|---|
| 340 | // printing all inserted for debugging
|
|---|
| 341 | // {
|
|---|
| 342 | // *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
|
|---|
| 343 | // int i=0;
|
|---|
| 344 | // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
|
|---|
| 345 | // if (runner != BoundaryPoints[axis].begin())
|
|---|
| 346 | // *out << ", " << i << ": " << *runner->second.second;
|
|---|
| 347 | // else
|
|---|
| 348 | // *out << i << ": " << *runner->second.second;
|
|---|
| 349 | // i++;
|
|---|
| 350 | // }
|
|---|
| 351 | // *out << endl;
|
|---|
| 352 | // }
|
|---|
| 353 | // 3c. throw out points whose distance is less than the mean of left and right neighbours
|
|---|
| 354 | bool flag = false;
|
|---|
| 355 | *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
|
|---|
| 356 | do { // do as long as we still throw one out per round
|
|---|
| 357 | flag = false;
|
|---|
| 358 | Boundaries::iterator left = BoundaryPoints[axis].end();
|
|---|
| 359 | Boundaries::iterator right = BoundaryPoints[axis].end();
|
|---|
| 360 | for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
|
|---|
| 361 | // set neighbours correctly
|
|---|
| 362 | if (runner == BoundaryPoints[axis].begin()) {
|
|---|
| 363 | left = BoundaryPoints[axis].end();
|
|---|
| 364 | } else {
|
|---|
| 365 | left = runner;
|
|---|
| 366 | }
|
|---|
| 367 | left--;
|
|---|
| 368 | right = runner;
|
|---|
| 369 | right++;
|
|---|
| 370 | if (right == BoundaryPoints[axis].end()) {
|
|---|
| 371 | right = BoundaryPoints[axis].begin();
|
|---|
| 372 | }
|
|---|
| 373 | // check distance
|
|---|
| 374 |
|
|---|
| 375 | // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
|
|---|
| 376 | {
|
|---|
| 377 | Vector SideA, SideB, SideC, SideH;
|
|---|
| 378 | SideA.CopyVector(&left->second.second->x);
|
|---|
| 379 | SideA.SubtractVector(MolCenter);
|
|---|
| 380 | SideA.ProjectOntoPlane(&AxisVector);
|
|---|
| 381 | // *out << "SideA: ";
|
|---|
| 382 | // SideA.Output(out);
|
|---|
| 383 | // *out << endl;
|
|---|
| 384 |
|
|---|
| 385 | SideB.CopyVector(&right->second.second->x);
|
|---|
| 386 | SideB.SubtractVector(MolCenter);
|
|---|
| 387 | SideB.ProjectOntoPlane(&AxisVector);
|
|---|
| 388 | // *out << "SideB: ";
|
|---|
| 389 | // SideB.Output(out);
|
|---|
| 390 | // *out << endl;
|
|---|
| 391 |
|
|---|
| 392 | SideC.CopyVector(&left->second.second->x);
|
|---|
| 393 | SideC.SubtractVector(&right->second.second->x);
|
|---|
| 394 | SideC.ProjectOntoPlane(&AxisVector);
|
|---|
| 395 | // *out << "SideC: ";
|
|---|
| 396 | // SideC.Output(out);
|
|---|
| 397 | // *out << endl;
|
|---|
| 398 |
|
|---|
| 399 | SideH.CopyVector(&runner->second.second->x);
|
|---|
| 400 | SideH.SubtractVector(MolCenter);
|
|---|
| 401 | SideH.ProjectOntoPlane(&AxisVector);
|
|---|
| 402 | // *out << "SideH: ";
|
|---|
| 403 | // SideH.Output(out);
|
|---|
| 404 | // *out << endl;
|
|---|
| 405 |
|
|---|
| 406 | // calculate each length
|
|---|
| 407 | double a = SideA.Norm();
|
|---|
| 408 | //double b = SideB.Norm();
|
|---|
| 409 | //double c = SideC.Norm();
|
|---|
| 410 | double h = SideH.Norm();
|
|---|
| 411 | // calculate the angles
|
|---|
| 412 | double alpha = SideA.Angle(&SideH);
|
|---|
| 413 | double beta = SideA.Angle(&SideC);
|
|---|
| 414 | double gamma = SideB.Angle(&SideH);
|
|---|
| 415 | double delta = SideC.Angle(&SideH);
|
|---|
| 416 | double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
|
|---|
| 417 | //*out << Verbose(2) << " 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;
|
|---|
| 418 | *out << 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;
|
|---|
| 419 | if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
|
|---|
| 420 | // throw out point
|
|---|
| 421 | *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
|
|---|
| 422 | BoundaryPoints[axis].erase(runner);
|
|---|
| 423 | flag = true;
|
|---|
| 424 | }
|
|---|
| 425 | }
|
|---|
| 426 | }
|
|---|
| 427 | } while (flag);
|
|---|
| 428 | }
|
|---|
| 429 | delete(MolCenter);
|
|---|
| 430 | return BoundaryPoints;
|
|---|
| 431 | };
|
|---|
| 432 |
|
|---|
| 433 | /** Tesselates the convex boundary by finding all boundary points.
|
|---|
| 434 | * \param *out output stream for debugging
|
|---|
| 435 | * \param *mol molecule structure with Atom's and Bond's
|
|---|
| 436 | * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
|
|---|
| 437 | * \param *LCList atoms in LinkedCell list
|
|---|
| 438 | * \param *filename filename prefix for output of vertex data
|
|---|
| 439 | * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
|
|---|
| 440 | */
|
|---|
| 441 | void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
|
|---|
| 442 | {
|
|---|
| 443 | bool BoundaryFreeFlag = false;
|
|---|
| 444 | Boundaries *BoundaryPoints = NULL;
|
|---|
| 445 |
|
|---|
| 446 | cout << Verbose(1) << "Begin of find_convex_border" << endl;
|
|---|
| 447 |
|
|---|
| 448 | if (TesselStruct != NULL) // free if allocated
|
|---|
| 449 | delete(TesselStruct);
|
|---|
| 450 | TesselStruct = new class Tesselation;
|
|---|
| 451 |
|
|---|
| 452 | // 1. Find all points on the boundary
|
|---|
| 453 | if (BoundaryPoints == NULL) {
|
|---|
| 454 | BoundaryFreeFlag = true;
|
|---|
| 455 | BoundaryPoints = GetBoundaryPoints(out, mol);
|
|---|
| 456 | } else {
|
|---|
| 457 | *out << Verbose(1) << "Using given boundary points set." << endl;
|
|---|
| 458 | }
|
|---|
| 459 |
|
|---|
| 460 | // printing all inserted for debugging
|
|---|
| 461 | for (int axis=0; axis < NDIM; axis++)
|
|---|
| 462 | {
|
|---|
| 463 | *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
|
|---|
| 464 | int i=0;
|
|---|
| 465 | for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
|
|---|
| 466 | if (runner != BoundaryPoints[axis].begin())
|
|---|
| 467 | *out << ", " << i << ": " << *runner->second.second;
|
|---|
| 468 | else
|
|---|
| 469 | *out << i << ": " << *runner->second.second;
|
|---|
| 470 | i++;
|
|---|
| 471 | }
|
|---|
| 472 | *out << endl;
|
|---|
| 473 | }
|
|---|
| 474 |
|
|---|
| 475 | // 2. fill the boundary point list
|
|---|
| 476 | for (int axis = 0; axis < NDIM; axis++)
|
|---|
| 477 | for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
|
|---|
| 478 | TesselStruct->AddPoint(runner->second.second);
|
|---|
| 479 |
|
|---|
| 480 | *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
|
|---|
| 481 | // now we have the whole set of edge points in the BoundaryList
|
|---|
| 482 |
|
|---|
| 483 | // listing for debugging
|
|---|
| 484 | // *out << Verbose(1) << "Listing PointsOnBoundary:";
|
|---|
| 485 | // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
|
|---|
| 486 | // *out << " " << *runner->second;
|
|---|
| 487 | // }
|
|---|
| 488 | // *out << endl;
|
|---|
| 489 |
|
|---|
| 490 | // 3a. guess starting triangle
|
|---|
| 491 | TesselStruct->GuessStartingTriangle(out);
|
|---|
| 492 |
|
|---|
| 493 | // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
|
|---|
| 494 | TesselStruct->TesselateOnBoundary(out, mol);
|
|---|
| 495 |
|
|---|
| 496 | // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
|
|---|
| 497 | if (!TesselStruct->InsertStraddlingPoints(out, mol))
|
|---|
| 498 | *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
|
|---|
| 499 |
|
|---|
| 500 | // 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
|
|---|
| 501 | if (!TesselStruct->CorrectConcaveBaselines(out))
|
|---|
| 502 | *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
|
|---|
| 503 |
|
|---|
| 504 | *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
|
|---|
| 505 |
|
|---|
| 506 | // 4. Store triangles in tecplot file
|
|---|
| 507 | if (filename != NULL) {
|
|---|
| 508 | if (DoTecplotOutput) {
|
|---|
| 509 | string OutputName(filename);
|
|---|
| 510 | OutputName.append(TecplotSuffix);
|
|---|
| 511 | ofstream *tecplot = new ofstream(OutputName.c_str());
|
|---|
| 512 | write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
|
|---|
| 513 | tecplot->close();
|
|---|
| 514 | delete(tecplot);
|
|---|
| 515 | }
|
|---|
| 516 | if (DoRaster3DOutput) {
|
|---|
| 517 | string OutputName(filename);
|
|---|
| 518 | OutputName.append(Raster3DSuffix);
|
|---|
| 519 | ofstream *rasterplot = new ofstream(OutputName.c_str());
|
|---|
| 520 | write_raster3d_file(out, rasterplot, TesselStruct, mol);
|
|---|
| 521 | rasterplot->close();
|
|---|
| 522 | delete(rasterplot);
|
|---|
| 523 | }
|
|---|
| 524 | }
|
|---|
| 525 |
|
|---|
| 526 | // free reference lists
|
|---|
| 527 | if (BoundaryFreeFlag)
|
|---|
| 528 | delete[] (BoundaryPoints);
|
|---|
| 529 |
|
|---|
| 530 | cout << Verbose(1) << "End of find_convex_border" << endl;
|
|---|
| 531 | };
|
|---|
| 532 |
|
|---|
| 533 |
|
|---|
| 534 | /** Determines the volume of a cluster.
|
|---|
| 535 | * Determines first the convex envelope, then tesselates it and calculates its volume.
|
|---|
| 536 | * \param *out output stream for debugging
|
|---|
| 537 | * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
|
|---|
| 538 | * \param *configuration needed for path to store convex envelope file
|
|---|
| 539 | * \return determined volume of the cluster in cubed config:GetIsAngstroem()
|
|---|
| 540 | */
|
|---|
| 541 | double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
|
|---|
| 542 | {
|
|---|
| 543 | bool IsAngstroem = configuration->GetIsAngstroem();
|
|---|
| 544 | double volume = 0.;
|
|---|
| 545 | double PyramidVolume = 0.;
|
|---|
| 546 | double G, h;
|
|---|
| 547 | Vector x, y;
|
|---|
| 548 | double a, b, c;
|
|---|
| 549 |
|
|---|
| 550 | // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
|
|---|
| 551 | *out << Verbose(1)
|
|---|
| 552 | << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
|
|---|
| 553 | << endl;
|
|---|
| 554 | for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
|
|---|
| 555 | { // go through every triangle, calculate volume of its pyramid with CoG as peak
|
|---|
| 556 | x.CopyVector(runner->second->endpoints[0]->node->node);
|
|---|
| 557 | x.SubtractVector(runner->second->endpoints[1]->node->node);
|
|---|
| 558 | y.CopyVector(runner->second->endpoints[0]->node->node);
|
|---|
| 559 | y.SubtractVector(runner->second->endpoints[2]->node->node);
|
|---|
| 560 | a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
|
|---|
| 561 | b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
|
|---|
| 562 | c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
|
|---|
| 563 | G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
|
|---|
| 564 | x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
|
|---|
| 565 | x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
|
|---|
| 566 | h = x.Norm(); // distance of CoG to triangle
|
|---|
| 567 | PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
|
|---|
| 568 | *out << Verbose(2) << "Area of triangle is " << G << " "
|
|---|
| 569 | << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
|
|---|
| 570 | << h << " and the volume is " << PyramidVolume << " "
|
|---|
| 571 | << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
|
|---|
| 572 | volume += PyramidVolume;
|
|---|
| 573 | }
|
|---|
| 574 | *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
|
|---|
| 575 | << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
|
|---|
| 576 | << endl;
|
|---|
| 577 |
|
|---|
| 578 | return volume;
|
|---|
| 579 | }
|
|---|
| 580 | ;
|
|---|
| 581 |
|
|---|
| 582 | /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
|
|---|
| 583 | * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
|
|---|
| 584 | * \param *out output stream for debugging
|
|---|
| 585 | * \param *configuration needed for path to store convex envelope file
|
|---|
| 586 | * \param *mol molecule structure representing the cluster
|
|---|
| 587 | * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
|
|---|
| 588 | * \param celldensity desired average density in final cell
|
|---|
| 589 | */
|
|---|
| 590 | void
|
|---|
| 591 | PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
|
|---|
| 592 | double ClusterVolume, double celldensity)
|
|---|
| 593 | {
|
|---|
| 594 | // transform to PAS
|
|---|
| 595 | mol->PrincipalAxisSystem(out, true);
|
|---|
| 596 |
|
|---|
| 597 | // some preparations beforehand
|
|---|
| 598 | bool IsAngstroem = configuration->GetIsAngstroem();
|
|---|
| 599 | Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
|
|---|
| 600 | class Tesselation *TesselStruct = NULL;
|
|---|
| 601 | LinkedCell LCList(mol, 10.);
|
|---|
| 602 | Find_convex_border(out, mol, TesselStruct, &LCList, NULL);
|
|---|
| 603 | double clustervolume;
|
|---|
| 604 | if (ClusterVolume == 0)
|
|---|
| 605 | clustervolume = VolumeOfConvexEnvelope(out, TesselStruct, configuration);
|
|---|
| 606 | else
|
|---|
| 607 | clustervolume = ClusterVolume;
|
|---|
| 608 | double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
|
|---|
| 609 | Vector BoxLengths;
|
|---|
| 610 | int repetition[NDIM] =
|
|---|
| 611 | { 1, 1, 1 };
|
|---|
| 612 | int TotalNoClusters = 1;
|
|---|
| 613 | for (int i = 0; i < NDIM; i++)
|
|---|
| 614 | TotalNoClusters *= repetition[i];
|
|---|
| 615 |
|
|---|
| 616 | // sum up the atomic masses
|
|---|
| 617 | double totalmass = 0.;
|
|---|
| 618 | atom *Walker = mol->start;
|
|---|
| 619 | while (Walker->next != mol->end)
|
|---|
| 620 | {
|
|---|
| 621 | Walker = Walker->next;
|
|---|
| 622 | totalmass += Walker->type->mass;
|
|---|
| 623 | }
|
|---|
| 624 | *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
|
|---|
| 625 | << totalmass << " atomicmassunit." << endl;
|
|---|
| 626 |
|
|---|
| 627 | *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
|
|---|
| 628 | << totalmass / clustervolume << " atomicmassunit/"
|
|---|
| 629 | << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
|
|---|
| 630 |
|
|---|
| 631 | // solve cubic polynomial
|
|---|
| 632 | *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
|
|---|
| 633 | << endl;
|
|---|
| 634 | double cellvolume;
|
|---|
| 635 | if (IsAngstroem)
|
|---|
| 636 | cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
|
|---|
| 637 | / clustervolume)) / (celldensity - 1);
|
|---|
| 638 | else
|
|---|
| 639 | cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
|
|---|
| 640 | / clustervolume)) / (celldensity - 1);
|
|---|
| 641 | *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
|
|---|
| 642 | << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
|
|---|
| 643 | : "atomiclength") << "^3." << endl;
|
|---|
| 644 |
|
|---|
| 645 | double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
|
|---|
| 646 | * GreatestDiameter[1] * GreatestDiameter[2]);
|
|---|
| 647 | *out << Verbose(1)
|
|---|
| 648 | << "Minimum volume of the convex envelope contained in a rectangular box is "
|
|---|
| 649 | << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
|
|---|
| 650 | : "atomiclength") << "^3." << endl;
|
|---|
| 651 | if (minimumvolume > cellvolume)
|
|---|
| 652 | {
|
|---|
| 653 | cerr << Verbose(0)
|
|---|
| 654 | << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
|
|---|
| 655 | << endl;
|
|---|
| 656 | cout << Verbose(0)
|
|---|
| 657 | << "Setting Box dimensions to minimum possible, the greatest diameters."
|
|---|
| 658 | << endl;
|
|---|
| 659 | for (int i = 0; i < NDIM; i++)
|
|---|
| 660 | BoxLengths.x[i] = GreatestDiameter[i];
|
|---|
| 661 | mol->CenterEdge(out, &BoxLengths);
|
|---|
| 662 | }
|
|---|
| 663 | else
|
|---|
| 664 | {
|
|---|
| 665 | BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
|
|---|
| 666 | * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
|
|---|
| 667 | BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
|
|---|
| 668 | * GreatestDiameter[1] + repetition[0] * repetition[2]
|
|---|
| 669 | * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
|
|---|
| 670 | * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
|
|---|
| 671 | BoxLengths.x[2] = minimumvolume - cellvolume;
|
|---|
| 672 | double x0 = 0., x1 = 0., x2 = 0.;
|
|---|
| 673 | if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
|
|---|
| 674 | BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
|
|---|
| 675 | *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
|
|---|
| 676 | << " ." << endl;
|
|---|
| 677 | else
|
|---|
| 678 | {
|
|---|
| 679 | *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
|
|---|
| 680 | << " and " << x1 << " and " << x2 << " ." << endl;
|
|---|
| 681 | x0 = x2; // sorted in ascending order
|
|---|
| 682 | }
|
|---|
| 683 |
|
|---|
| 684 | cellvolume = 1;
|
|---|
| 685 | for (int i = 0; i < NDIM; i++)
|
|---|
| 686 | {
|
|---|
| 687 | BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
|
|---|
| 688 | cellvolume *= BoxLengths.x[i];
|
|---|
| 689 | }
|
|---|
| 690 |
|
|---|
| 691 | // set new box dimensions
|
|---|
| 692 | *out << Verbose(0) << "Translating to box with these boundaries." << endl;
|
|---|
| 693 | mol->SetBoxDimension(&BoxLengths);
|
|---|
| 694 | mol->CenterInBox((ofstream *) &cout);
|
|---|
| 695 | }
|
|---|
| 696 | // update Box of atoms by boundary
|
|---|
| 697 | mol->SetBoxDimension(&BoxLengths);
|
|---|
| 698 | *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
|
|---|
| 699 | << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
|
|---|
| 700 | << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
|
|---|
| 701 | << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
|
|---|
| 702 | }
|
|---|
| 703 | ;
|
|---|
| 704 |
|
|---|
| 705 |
|
|---|
| 706 | /** Fills the empty space of the simulation box with water/
|
|---|
| 707 | * \param *out output stream for debugging
|
|---|
| 708 | * \param *List list of molecules already present in box
|
|---|
| 709 | * \param *TesselStruct contains tesselated surface
|
|---|
| 710 | * \param *filler molecule which the box is to be filled with
|
|---|
| 711 | * \param configuration contains box dimensions
|
|---|
| 712 | * \param distance[NDIM] distance between filling molecules in each direction
|
|---|
| 713 | * \param RandAtomDisplacement maximum distance for random displacement per atom
|
|---|
| 714 | * \param RandMolDisplacement maximum distance for random displacement per filler molecule
|
|---|
| 715 | * \param DoRandomRotation true - do random rotiations, false - don't
|
|---|
| 716 | * \return *mol pointer to new molecule with filled atoms
|
|---|
| 717 | */
|
|---|
| 718 | molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
|
|---|
| 719 | {
|
|---|
| 720 | molecule *Filling = new molecule(filler->elemente);
|
|---|
| 721 | Vector CurrentPosition;
|
|---|
| 722 | int N[NDIM];
|
|---|
| 723 | int n[NDIM];
|
|---|
| 724 | double *M = filler->ReturnFullMatrixforSymmetric(filler->cell_size);
|
|---|
| 725 | double Rotations[NDIM*NDIM];
|
|---|
| 726 | Vector AtomTranslations;
|
|---|
| 727 | Vector FillerTranslations;
|
|---|
| 728 | Vector FillerDistance;
|
|---|
| 729 | double FillIt = false;
|
|---|
| 730 | atom *Walker = NULL, *Runner = NULL;
|
|---|
| 731 | bond *Binder = NULL;
|
|---|
| 732 |
|
|---|
| 733 | // Center filler at origin
|
|---|
| 734 | filler->CenterOrigin(out);
|
|---|
| 735 | filler->Center.Zero();
|
|---|
| 736 |
|
|---|
| 737 | // calculate filler grid in [0,1]^3
|
|---|
| 738 | FillerDistance.Init(distance[0], distance[1], distance[2]);
|
|---|
| 739 | FillerDistance.InverseMatrixMultiplication(M);
|
|---|
| 740 | for(int i=0;i<NDIM;i++)
|
|---|
| 741 | N[i] = (int) ceil(1./FillerDistance.x[i]);
|
|---|
| 742 |
|
|---|
| 743 | // go over [0,1]^3 filler grid
|
|---|
| 744 | for (n[0] = 0; n[0] < N[0]; n[0]++)
|
|---|
| 745 | for (n[1] = 0; n[1] < N[1]; n[1]++)
|
|---|
| 746 | for (n[2] = 0; n[2] < N[2]; n[2]++) {
|
|---|
| 747 | // calculate position of current grid vector in untransformed box
|
|---|
| 748 | CurrentPosition.Init(n[0], n[1], n[2]);
|
|---|
| 749 | CurrentPosition.MatrixMultiplication(M);
|
|---|
| 750 | // Check whether point is in- or outside
|
|---|
| 751 | FillIt = true;
|
|---|
| 752 | for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
|
|---|
| 753 | FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
|
|---|
| 754 | }
|
|---|
| 755 |
|
|---|
| 756 | if (FillIt) {
|
|---|
| 757 | // fill in Filler
|
|---|
| 758 | *out << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
|
|---|
| 759 |
|
|---|
| 760 | // create molecule random translation vector ...
|
|---|
| 761 | for (int i=0;i<NDIM;i++)
|
|---|
| 762 | FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
|
|---|
| 763 |
|
|---|
| 764 | // go through all atoms
|
|---|
| 765 | Walker = filler->start;
|
|---|
| 766 | while (Walker != filler->end) {
|
|---|
| 767 | Walker = Walker->next;
|
|---|
| 768 | // copy atom ...
|
|---|
| 769 | *Runner = new atom(Walker);
|
|---|
| 770 |
|
|---|
| 771 | // create atomic random translation vector ...
|
|---|
| 772 | for (int i=0;i<NDIM;i++)
|
|---|
| 773 | AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
|
|---|
| 774 |
|
|---|
| 775 | // ... and rotation matrix
|
|---|
| 776 | if (DoRandomRotation) {
|
|---|
| 777 | double phi[NDIM];
|
|---|
| 778 | for (int i=0;i<NDIM;i++) {
|
|---|
| 779 | phi[i] = rand()/(RAND_MAX/(2.*M_PI));
|
|---|
| 780 | }
|
|---|
| 781 |
|
|---|
| 782 | Rotations[0] = cos(phi[0]) *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
|
|---|
| 783 | Rotations[3] = sin(phi[0]) *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
|
|---|
| 784 | Rotations[6] = cos(phi[1])*sin(phi[2]) ;
|
|---|
| 785 | Rotations[1] = - sin(phi[0])*cos(phi[1]) ;
|
|---|
| 786 | Rotations[4] = cos(phi[0])*cos(phi[1]) ;
|
|---|
| 787 | Rotations[7] = sin(phi[1]) ;
|
|---|
| 788 | Rotations[3] = - cos(phi[0]) *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
|
|---|
| 789 | Rotations[5] = - sin(phi[0]) *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
|
|---|
| 790 | Rotations[8] = cos(phi[1])*cos(phi[2]) ;
|
|---|
| 791 | }
|
|---|
| 792 |
|
|---|
| 793 | // ... and put at new position
|
|---|
| 794 | if (DoRandomRotation)
|
|---|
| 795 | Runner->x.MatrixMultiplication(Rotations);
|
|---|
| 796 | Runner->x.AddVector(&AtomTranslations);
|
|---|
| 797 | Runner->x.AddVector(&FillerTranslations);
|
|---|
| 798 | Runner->x.AddVector(&CurrentPosition);
|
|---|
| 799 | // insert into Filling
|
|---|
| 800 | Filling->AddAtom(Runner);
|
|---|
| 801 | }
|
|---|
| 802 |
|
|---|
| 803 | // go through all bonds and add as well
|
|---|
| 804 | Binder = filler->first;
|
|---|
| 805 | while(Binder != filler->last) {
|
|---|
| 806 | Binder = Binder->next;
|
|---|
| 807 | }
|
|---|
| 808 | } else {
|
|---|
| 809 | // leave empty
|
|---|
| 810 | *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
|
|---|
| 811 | }
|
|---|
| 812 | }
|
|---|
| 813 | return Filling;
|
|---|
| 814 | };
|
|---|
| 815 |
|
|---|
| 816 |
|
|---|
| 817 | /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
|
|---|
| 818 | * \param *out output stream for debugging
|
|---|
| 819 | * \param *mol molecule structure with Atom's and Bond's
|
|---|
| 820 | * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
|
|---|
| 821 | * \param *LCList atoms in LinkedCell list
|
|---|
| 822 | * \param *filename filename prefix for output of vertex data
|
|---|
| 823 | * \para RADIUS radius of the virtual sphere
|
|---|
| 824 | */
|
|---|
| 825 | void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
|
|---|
| 826 | {
|
|---|
| 827 | int N = 0;
|
|---|
| 828 | bool freeTess = false;
|
|---|
| 829 | bool freeLC = false;
|
|---|
| 830 | ofstream *tempstream = NULL;
|
|---|
| 831 | char NumberName[255];
|
|---|
| 832 | int TriangleFilesWritten = 0;
|
|---|
| 833 |
|
|---|
| 834 | *out << Verbose(1) << "Entering search for non convex hull. " << endl;
|
|---|
| 835 | if (Tess == NULL) {
|
|---|
| 836 | *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
|
|---|
| 837 | Tess = new Tesselation;
|
|---|
| 838 | freeTess = true;
|
|---|
| 839 | }
|
|---|
| 840 | LineMap::iterator baseline;
|
|---|
| 841 | LineMap::iterator testline;
|
|---|
| 842 | *out << Verbose(0) << "Begin of Find_non_convex_border\n";
|
|---|
| 843 | bool flag = false; // marks whether we went once through all baselines without finding any without two triangles
|
|---|
| 844 | bool failflag = false;
|
|---|
| 845 |
|
|---|
| 846 | if (LCList == NULL) {
|
|---|
| 847 | LCList = new LinkedCell(mol, 2.*RADIUS);
|
|---|
| 848 | freeLC = true;
|
|---|
| 849 | }
|
|---|
| 850 |
|
|---|
| 851 | Tess->Find_starting_triangle(out, RADIUS, LCList);
|
|---|
| 852 |
|
|---|
| 853 | baseline = Tess->LinesOnBoundary.begin();
|
|---|
| 854 | while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
|
|---|
| 855 | if (baseline->second->TrianglesCount == 1) {
|
|---|
| 856 | failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
|
|---|
| 857 | flag = flag || failflag;
|
|---|
| 858 | if (!failflag)
|
|---|
| 859 | cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
|
|---|
| 860 | } else {
|
|---|
| 861 | //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
|
|---|
| 862 | if (baseline->second->TrianglesCount != 2)
|
|---|
| 863 | cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
|
|---|
| 864 | }
|
|---|
| 865 |
|
|---|
| 866 | N++;
|
|---|
| 867 | baseline++;
|
|---|
| 868 | if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
|
|---|
| 869 | baseline = Tess->LinesOnBoundary.begin(); // restart if we reach end due to newly inserted lines
|
|---|
| 870 | flag = false;
|
|---|
| 871 | }
|
|---|
| 872 |
|
|---|
| 873 | // write temporary envelope
|
|---|
| 874 | if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
|
|---|
| 875 | class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
|
|---|
| 876 | sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
|
|---|
| 877 | if (DoTecplotOutput) {
|
|---|
| 878 | string NameofTempFile(filename);
|
|---|
| 879 | NameofTempFile.append(NumberName);
|
|---|
| 880 | for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
|
|---|
| 881 | NameofTempFile.erase(npos, 1);
|
|---|
| 882 | NameofTempFile.append(TecplotSuffix);
|
|---|
| 883 | cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
|
|---|
| 884 | tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
|
|---|
| 885 | write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
|
|---|
| 886 | tempstream->close();
|
|---|
| 887 | tempstream->flush();
|
|---|
| 888 | delete(tempstream);
|
|---|
| 889 | }
|
|---|
| 890 |
|
|---|
| 891 | if (DoRaster3DOutput) {
|
|---|
| 892 | string NameofTempFile(filename);
|
|---|
| 893 | NameofTempFile.append(NumberName);
|
|---|
| 894 | for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
|
|---|
| 895 | NameofTempFile.erase(npos, 1);
|
|---|
| 896 | NameofTempFile.append(Raster3DSuffix);
|
|---|
| 897 | cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
|
|---|
| 898 | tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
|
|---|
| 899 | write_raster3d_file(out, tempstream, Tess, mol);
|
|---|
| 900 | // // include the current position of the virtual sphere in the temporary raster3d file
|
|---|
| 901 | // // make the circumsphere's center absolute again
|
|---|
| 902 | // helper.CopyVector(BaseRay->endpoints[0]->node->node);
|
|---|
| 903 | // helper.AddVector(BaseRay->endpoints[1]->node->node);
|
|---|
| 904 | // helper.Scale(0.5);
|
|---|
| 905 | // (*it)->OptCenter.AddVector(&helper);
|
|---|
| 906 | // Vector *center = mol->DetermineCenterOfAll(out);
|
|---|
| 907 | // (*it)->OptCenter.SubtractVector(center);
|
|---|
| 908 | // delete(center);
|
|---|
| 909 | // // and add to file plus translucency object
|
|---|
| 910 | // *tempstream << "# current virtual sphere\n";
|
|---|
| 911 | // *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n";
|
|---|
| 912 | // *tempstream << "2\n " << (*it)->OptCenter.x[0] << " "
|
|---|
| 913 | // << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
|
|---|
| 914 | // << "\t" << RADIUS << "\t1 0 0\n";
|
|---|
| 915 | // *tempstream << "9\n terminating special property\n";
|
|---|
| 916 | tempstream->close();
|
|---|
| 917 | tempstream->flush();
|
|---|
| 918 | delete(tempstream);
|
|---|
| 919 | }
|
|---|
| 920 | if (DoTecplotOutput || DoRaster3DOutput)
|
|---|
| 921 | TriangleFilesWritten++;
|
|---|
| 922 | }
|
|---|
| 923 | }
|
|---|
| 924 |
|
|---|
| 925 | // write final envelope
|
|---|
| 926 | if (filename != 0) {
|
|---|
| 927 | *out << Verbose(1) << "Writing final tecplot file\n";
|
|---|
| 928 | if (DoTecplotOutput) {
|
|---|
| 929 | string OutputName(filename);
|
|---|
| 930 | OutputName.append(TecplotSuffix);
|
|---|
| 931 | ofstream *tecplot = new ofstream(OutputName.c_str());
|
|---|
| 932 | write_tecplot_file(out, tecplot, Tess, mol, -1);
|
|---|
| 933 | tecplot->close();
|
|---|
| 934 | delete(tecplot);
|
|---|
| 935 | }
|
|---|
| 936 | if (DoRaster3DOutput) {
|
|---|
| 937 | string OutputName(filename);
|
|---|
| 938 | OutputName.append(Raster3DSuffix);
|
|---|
| 939 | ofstream *raster = new ofstream(OutputName.c_str());
|
|---|
| 940 | write_raster3d_file(out, raster, Tess, mol);
|
|---|
| 941 | raster->close();
|
|---|
| 942 | delete(raster);
|
|---|
| 943 | }
|
|---|
| 944 | } else {
|
|---|
| 945 | cerr << "ERROR: Could definitively not find all necessary triangles!" << endl;
|
|---|
| 946 | }
|
|---|
| 947 |
|
|---|
| 948 | cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
|
|---|
| 949 | int counter = 0;
|
|---|
| 950 | for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) {
|
|---|
| 951 | if (testline->second->TrianglesCount != 2) {
|
|---|
| 952 | cout << Verbose(2) << *testline->second << "\t" << testline->second->TrianglesCount << endl;
|
|---|
| 953 | counter++;
|
|---|
| 954 | }
|
|---|
| 955 | }
|
|---|
| 956 | if (counter == 0)
|
|---|
| 957 | cout << Verbose(2) << "None." << endl;
|
|---|
| 958 |
|
|---|
| 959 | // Tests the IsInnerAtom() function.
|
|---|
| 960 | Vector x (0, 0, 0);
|
|---|
| 961 | cout << Verbose(0) << "Point to check: " << x << endl;
|
|---|
| 962 | cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
|
|---|
| 963 | << "for vector " << x << "." << endl;
|
|---|
| 964 | TesselPoint* a = Tess->PointsOnBoundary.begin()->second->node;
|
|---|
| 965 | cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
|
|---|
| 966 | cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerPoint(a, Tess, LCList)
|
|---|
| 967 | << "for atom " << a << " (on boundary)." << endl;
|
|---|
| 968 | LinkedNodes *List = NULL;
|
|---|
| 969 | for (int i=0;i<NDIM;i++) { // each axis
|
|---|
| 970 | LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
|
|---|
| 971 | for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
|
|---|
| 972 | for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
|
|---|
| 973 | List = LCList->GetCurrentCell();
|
|---|
| 974 | //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
|
|---|
| 975 | if (List != NULL) {
|
|---|
| 976 | for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
|
|---|
| 977 | if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
|
|---|
| 978 | a = *Runner;
|
|---|
| 979 | i=3;
|
|---|
| 980 | for (int j=0;j<NDIM;j++)
|
|---|
| 981 | LCList->n[j] = LCList->N[j];
|
|---|
| 982 | break;
|
|---|
| 983 | }
|
|---|
| 984 | }
|
|---|
| 985 | }
|
|---|
| 986 | }
|
|---|
| 987 | }
|
|---|
| 988 | cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(a, Tess, LCList)
|
|---|
| 989 | << "for atom " << a << " (inside)." << endl;
|
|---|
| 990 |
|
|---|
| 991 |
|
|---|
| 992 | if (freeTess)
|
|---|
| 993 | delete(Tess);
|
|---|
| 994 | if (freeLC)
|
|---|
| 995 | delete(LCList);
|
|---|
| 996 | *out << Verbose(0) << "End of Find_non_convex_border\n";
|
|---|
| 997 | };
|
|---|
| 998 |
|
|---|
| 999 | /** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
|
|---|
| 1000 | * \param *out output stream for debugging
|
|---|
| 1001 | * \param *srcmol molecule to embed into
|
|---|
| 1002 | * \return *Vector new center of \a *srcmol for embedding relative to \a this
|
|---|
| 1003 | */
|
|---|
| 1004 | Vector* molecule::FindEmbeddingHole(ofstream *out, molecule *srcmol)
|
|---|
| 1005 | {
|
|---|
| 1006 | Vector *Center = new Vector;
|
|---|
| 1007 | Center->Zero();
|
|---|
| 1008 | // calculate volume/shape of \a *srcmol
|
|---|
| 1009 |
|
|---|
| 1010 | // find embedding holes
|
|---|
| 1011 |
|
|---|
| 1012 | // if more than one, let user choose
|
|---|
| 1013 |
|
|---|
| 1014 | // return embedding center
|
|---|
| 1015 | return Center;
|
|---|
| 1016 | };
|
|---|
| 1017 |
|
|---|
