Changes in / [51c910:36ec71]
- Location:
- src
- Files:
-
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundary.cpp
r51c910 r36ec71 1 1 #include "boundary.hpp" 2 #include "linkedcell.hpp" 3 #include "molecules.hpp" 4 #include <gsl/gsl_matrix.h> 5 #include <gsl/gsl_linalg.h> 6 #include <gsl/gsl_multimin.h> 7 #include <gsl/gsl_permutation.h> 2 8 3 9 #define DEBUG 1 … … 9 15 #define Raster3DSuffix ".r3d" 10 16 #define VRMLSUffix ".wrl" 11 #define HULLEPSILON MYEPSILON17 #define HULLEPSILON 1e-7 12 18 13 19 // ======================================== Points on Boundary ================================= … … 34 40 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl; 35 41 node = NULL; 36 lines.clear();37 42 } 38 43 ; … … 90 95 BoundaryLineSet::~BoundaryLineSet() 91 96 { 97 int Numbers[2]; 98 Numbers[0] = endpoints[1]->Nr; 99 Numbers[1] = endpoints[0]->Nr; 92 100 for (int i = 0; i < 2; i++) { 93 101 cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 94 endpoints[i]->lines.erase(Nr); 95 LineMap::iterator tester = endpoints[i]->lines.begin(); 96 tester++; 97 if (tester == endpoints[i]->lines.end()) { 102 // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set 103 pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]); 104 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 105 if ((*Runner).second == this) { 106 endpoints[i]->lines.erase(Runner); 107 break; 108 } 109 if (endpoints[i]->lines.empty()) { 98 110 cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl; 99 111 if (endpoints[i] != NULL) { … … 143 155 ; 144 156 145 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], 146 int number) 157 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) 147 158 { 148 159 // set number … … 191 202 lines[i]->triangles.erase(Nr); 192 203 if (lines[i]->triangles.empty()) { 193 cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;194 204 if (lines[i] != NULL) { 205 cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl; 195 206 delete (lines[i]); 196 207 lines[i] = NULL; … … 198 209 cerr << "ERROR: This line " << i << " has already been free'd." << endl; 199 210 } else 200 cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;211 cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl; 201 212 } 202 213 } … … 224 235 } 225 236 ; 237 238 239 // ============================ CandidateForTesselation ============================= 240 241 CandidateForTesselation::CandidateForTesselation( 242 atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter 243 ) { 244 point = candidate; 245 BaseLine = line; 246 OptCenter.CopyVector(&OptCandidateCenter); 247 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); 248 } 249 250 CandidateForTesselation::~CandidateForTesselation() { 251 point = NULL; 252 BaseLine = NULL; 253 } 226 254 227 255 // ========================================== F U N C T I O N S ================================= … … 672 700 } 673 701 *rasterfile << "1. 0. 0." << endl; // red as colour 674 *rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object702 //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object 675 703 } 676 704 *rasterfile << "9\n terminating special property\n"; … … 983 1011 // set new box dimensions 984 1012 *out << Verbose(0) << "Translating to box with these boundaries." << endl; 985 mol->CenterInBox((ofstream *) &cout, &BoxLengths); 1013 mol->SetBoxDimension(&BoxLengths); 1014 mol->CenterInBox((ofstream *) &cout); 986 1015 } 987 1016 // update Box of atoms by boundary … … 1019 1048 } else 1020 1049 cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl; 1021 }1022 for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {1023 if (runner->second != NULL) {1024 delete (runner->second);1025 runner->second = NULL;1026 } else1027 cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;1028 }1029 for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {1030 if (runner->second != NULL) {1031 delete (runner->second);1032 runner->second = NULL;1033 } else1034 cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;1035 1050 } 1036 1051 } … … 1464 1479 TPS[n] = new class BoundaryPointSet(Candidate); 1465 1480 InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n])); 1466 if (InsertUnique.second) // if new point was not present before, increase counter 1467 { 1468 PointsOnBoundaryCount++; 1469 } 1470 else 1471 { 1472 delete TPS[n]; 1473 cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node) 1474 << " gibt's schon in der PointMap." << endl; 1475 TPS[n] = (InsertUnique.first)->second; 1476 } 1481 if (InsertUnique.second) { // if new point was not present before, increase counter 1482 PointsOnBoundaryCount++; 1483 } else { 1484 delete TPS[n]; 1485 cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl; 1486 TPS[n] = (InsertUnique.first)->second; 1487 } 1477 1488 } 1478 1489 ; 1479 1490 1480 1491 /** Function tries to add line from current Points in BPS to BoundaryLineSet. 1481 * If succes ful it raises the line count and inserts the new line into the BLS,1482 * if unsucces ful, it writes the line which had been present into the BLS, deleting the new constructed one.1492 * If successful it raises the line count and inserts the new line into the BLS, 1493 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one. 1483 1494 * @param *a first endpoint 1484 1495 * @param *b second endpoint 1485 1496 * @param n index of Tesselation::BLS giving the line with both endpoints 1486 1497 */ 1487 void 1488 Tesselation::AddTriangleLine(class BoundaryPointSet *a, 1489 class BoundaryPointSet *b, int n) 1490 { 1491 LineMap::iterator LineWalker; 1492 //cout << "Manually checking endpoints for line." << endl; 1493 if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr) 1494 //If a line is there, how do I recognize that beyond a shadow of a doubt? 1495 { 1496 //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl; 1497 1498 LineWalker = LinesOnBoundary.end(); 1499 LineWalker--; 1500 1501 while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr, 1502 b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max( 1503 a->node->nr, b->node->nr)) 1504 { 1505 //cout << Verbose(1) << "Looking for line which already exists"<< endl; 1506 LineWalker--; 1507 } 1508 BPS[0] = LineWalker->second->endpoints[0]; 1509 BPS[1] = LineWalker->second->endpoints[1]; 1510 BLS[n] = LineWalker->second; 1511 1512 } 1513 else 1514 { 1515 cout << Verbose(2) 1516 << "Adding line which has not been used before between " 1517 << *(a->node) << " and " << *(b->node) << "." << endl; 1518 BPS[0] = a; 1519 BPS[1] = b; 1520 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1521 1522 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n])); 1523 LinesOnBoundaryCount++; 1524 1525 } 1498 void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) { 1499 bool insertNewLine = true; 1500 1501 if (a->lines.find(b->node->nr) != a->lines.end()) { 1502 LineMap::iterator FindLine; 1503 pair<LineMap::iterator,LineMap::iterator> FindPair; 1504 FindPair = a->lines.equal_range(b->node->nr); 1505 1506 for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) { 1507 // If there is a line with less than two attached triangles, we don't need a new line. 1508 if (FindLine->second->TrianglesCount < 2) { 1509 insertNewLine = false; 1510 cout << Verbose(3) << "Using existing line " << *FindLine->second << endl; 1511 1512 BPS[0] = FindLine->second->endpoints[0]; 1513 BPS[1] = FindLine->second->endpoints[1]; 1514 BLS[n] = FindLine->second; 1515 1516 break; 1517 } 1518 } 1519 } 1520 1521 if (insertNewLine) { 1522 AlwaysAddTriangleLine(a, b, n); 1523 } 1526 1524 } 1527 1525 ; 1526 1527 /** 1528 * Adds lines from each of the current points in the BPS to BoundaryLineSet. 1529 * Raises the line count and inserts the new line into the BLS. 1530 * 1531 * @param *a first endpoint 1532 * @param *b second endpoint 1533 * @param n index of Tesselation::BLS giving the line with both endpoints 1534 */ 1535 void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) 1536 { 1537 cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl; 1538 BPS[0] = a; 1539 BPS[1] = b; 1540 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps 1541 // add line to global map 1542 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n])); 1543 // increase counter 1544 LinesOnBoundaryCount++; 1545 }; 1528 1546 1529 1547 /** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm). … … 1531 1549 */ 1532 1550 void 1533 Tesselation::AddTriangleToLines() 1534 { 1535 1536 cout << Verbose(1) << "Adding triangle to its lines" << endl; 1551 Tesselation::AddTriangle() 1552 { 1553 cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1554 1555 // add triangle to global map 1537 1556 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1538 1557 TrianglesOnBoundaryCount++; 1539 1558 1540 /* 1541 * this is apparently done when constructing triangle 1542 1543 for (int i=0; i<3; i++) 1544 { 1545 BLS[i]->AddTriangle(BTS); 1546 } 1547 */ 1559 // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet 1548 1560 } 1549 1561 ; … … 1632 1644 * @param b vector second point of triangle 1633 1645 * @param c vector third point of triangle 1646 * @param *Umkreismittelpunkt new cneter point of circumference 1634 1647 * @param Direction vector indicates up/down 1635 1648 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle … … 1685 1698 } 1686 1699 } 1700 1687 1701 TempNormal.Normalize(); 1688 1702 Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius); … … 1697 1711 cout << Verbose(3) << "End of Get_center_of_sphere.\n"; 1698 1712 }; 1699 1700 /** This recursive function finds a third point, to form a triangle with two given ones.1701 * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \1702 * supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \1703 * upon which we operate.1704 * If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \1705 * direction and angle into Storage.1706 * We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \1707 * with all neighbours of the candidate.1708 * @param a first point1709 * @param b second point1710 * *param c atom old third point of old triangle1711 * @param Candidate base point along whose bonds to start looking from1712 * @param Parent point to avoid during search as its wrong direction1713 * @param RecursionLevel contains current recursion depth1714 * @param Chord baseline vector of first and second point1715 * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to1716 * @param OldNormal normal of the triangle which the baseline belongs to1717 * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere1718 * @param Opt_Candidate candidate reference to return1719 * @param Storage array containing two angles of current Opt_Candidate1720 * @param RADIUS radius of ball1721 * @param mol molecule structure with atoms and bonds1722 */1723 void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,1724 int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,1725 atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)1726 {1727 cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";1728 cout << Verbose(3) << "Candidate is "<< *Candidate << endl;1729 cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;1730 cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;1731 cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;1732 cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;1733 /* OldNormal is normal vector on the old triangle1734 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.1735 * Chord points from b to a!!!1736 */1737 Vector dif_a; //Vector from a to candidate1738 Vector dif_b; //Vector from b to candidate1739 Vector AngleCheck;1740 Vector TempNormal, Umkreismittelpunkt;1741 Vector Mittelpunkt;1742 1743 double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;1744 double BallAngle, AlternativeSign;1745 atom *Walker; // variable atom point1746 1747 Vector NewUmkreismittelpunkt;1748 1749 if (a != Candidate and b != Candidate and c != Candidate) {1750 cout << Verbose(3) << "We have a unique candidate!" << endl;1751 dif_a.CopyVector(&(a->x));1752 dif_a.SubtractVector(&(Candidate->x));1753 dif_b.CopyVector(&(b->x));1754 dif_b.SubtractVector(&(Candidate->x));1755 AngleCheck.CopyVector(&(Candidate->x));1756 AngleCheck.SubtractVector(&(a->x));1757 AngleCheck.ProjectOntoPlane(Chord);1758 1759 SideA = dif_b.Norm();1760 SideB = dif_a.Norm();1761 SideC = Chord->Norm();1762 //Chord->Scale(-1);1763 1764 alpha = Chord->Angle(&dif_a);1765 beta = M_PI - Chord->Angle(&dif_b);1766 gamma = dif_a.Angle(&dif_b);1767 1768 cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;1769 1770 if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {1771 cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";1772 cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;1773 }1774 1775 if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {1776 Umkreisradius = SideA / 2.0 / sin(alpha);1777 //cout << Umkreisradius << endl;1778 //cout << SideB / 2.0 / sin(beta) << endl;1779 //cout << SideC / 2.0 / sin(gamma) << endl;1780 1781 if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.1782 cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;1783 cout << Verbose(2) << "Candidate is "<< *Candidate << endl;1784 sign = AngleCheck.ScalarProduct(direction1);1785 if (fabs(sign)<MYEPSILON) {1786 if (AngleCheck.ScalarProduct(OldNormal)<0) {1787 sign =0;1788 AlternativeSign=1;1789 } else {1790 sign =0;1791 AlternativeSign=-1;1792 }1793 } else {1794 sign /= fabs(sign);1795 cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;1796 }1797 if (sign >= 0) {1798 cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;1799 Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);1800 Mittelpunkt.SubtractVector(&ReferencePoint);1801 cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;1802 BallAngle = Mittelpunkt.Angle(OldNormal);1803 cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;1804 1805 1806 cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;1807 1808 if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {1809 if (Storage[0]< -1.5) { // first Candidate at all1810 if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {1811 cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";1812 Opt_Candidate = Candidate;1813 Storage[0] = sign;1814 Storage[1] = AlternativeSign;1815 Storage[2] = BallAngle;1816 cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;1817 } else1818 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;1819 } else {1820 if ( Storage[2] > BallAngle) {1821 if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {1822 cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";1823 Opt_Candidate = Candidate;1824 Storage[0] = sign;1825 Storage[1] = AlternativeSign;1826 Storage[2] = BallAngle;1827 cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;1828 } else1829 cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;1830 } else {1831 if (DEBUG) {1832 cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;1833 }1834 }1835 }1836 } else {1837 if (DEBUG) {1838 cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;1839 }1840 }1841 } else {1842 if (DEBUG) {1843 cout << Verbose(3) << *Candidate << " is not in search direction." << endl;1844 }1845 }1846 } else {1847 if (DEBUG) {1848 cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;1849 }1850 }1851 } else {1852 if (DEBUG) {1853 cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;1854 }1855 }1856 } else {1857 if (DEBUG) {1858 cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;1859 }1860 }1861 1862 if (RecursionLevel < 5) { // Seven is the recursion level threshold.1863 for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond1864 Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);1865 if (Walker == Parent) { // don't go back the same bond1866 continue;1867 } else {1868 Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again1869 }1870 }1871 }1872 cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";1873 }1874 ;1875 1713 1876 1714 … … 1895 1733 beta = M_PI - SideC.Angle(&SideA); 1896 1734 gamma = M_PI - SideA.Angle(&SideB); 1897 cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;1735 //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl; 1898 1736 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) 1899 cerr << " Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;1737 cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl; 1900 1738 1901 1739 Center->Zero(); … … 1944 1782 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals 1945 1783 alpha = 2.*M_PI - alpha; 1946 cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;1784 //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl; 1947 1785 radius = helper.Distance(&OldSphereCenter); 1948 1786 helper.ProjectOntoPlane(&NormalVector); 1949 1787 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles 1950 1788 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) { 1951 cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;1789 //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl; 1952 1790 return alpha; 1953 1791 } else { 1954 cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;1792 //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl; 1955 1793 return 2.*M_PI; 1956 1794 } … … 1958 1796 1959 1797 1960 /** This recursive function finds a third point, to form a triangle with two given ones. 1961 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points 1962 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then 1963 * the center of the sphere is still fixed up to a single parameter. The band of possible values 1964 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives 1965 * us the "null" on this circle, the new center of the candidate point will be some way along this 1966 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given 1967 * by the normal vector of the base triangle that always points outwards by construction. 1968 * Hence, we construct a Center of this circle which sits right in the middle of the current base line. 1969 * We construct the normal vector that defines the plane this circle lies in, it is just in the 1970 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest 1971 * with respect to the length of the baseline and the sphere's fixed \a RADIUS. 1972 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center 1973 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check 1974 * both. 1975 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check 1976 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check 1977 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal 1978 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality 1979 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether 1980 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI). 1981 * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points 1982 * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line 1983 * @param OptCandidate candidate reference on return 1984 * @param OptCandidateCenter candidate's sphere center on return 1985 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate 1986 * @param RADIUS radius of sphere 1987 * @param *LC LinkedCell structure with neighbouring atoms 1988 */ 1989 // void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC) 1990 // { 1991 // atom *Walker = NULL; 1992 // Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 1993 // Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 1994 // Vector OldSphereCenter; // center of the sphere defined by the three points of BaseTriangle 1995 // Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 1996 // Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 1997 // Vector NewNormalVector; // normal vector of the Candidate's triangle 1998 // Vector SearchDirection; // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate) 1999 // Vector helper; 2000 // LinkedAtoms *List = NULL; 2001 // double CircleRadius; // radius of this circle 2002 // double radius; 2003 // double alpha, Otheralpha; // angles (i.e. parameter for the circle). 2004 // double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle 2005 // int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2006 // atom *Candidate = NULL; 2007 // 2008 // cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl; 2009 // 2010 // cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl; 2011 // 2012 // // construct center of circle 2013 // CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x)); 2014 // CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x); 2015 // CircleCenter.Scale(0.5); 2016 // 2017 // // construct normal vector of circle 2018 // CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x); 2019 // CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x); 2020 // 2021 // // calculate squared radius of circle 2022 // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2023 // if (radius/4. < RADIUS*RADIUS) { 2024 // CircleRadius = RADIUS*RADIUS - radius/4.; 2025 // CirclePlaneNormal.Normalize(); 2026 // cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2027 // 2028 // // construct old center 2029 // GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x)); 2030 // helper.CopyVector(&BaseTriangle->NormalVector); // normal vector ensures that this is correct center of the two possible ones 2031 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter); 2032 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2033 // OldSphereCenter.AddVector(&helper); 2034 // OldSphereCenter.SubtractVector(&CircleCenter); 2035 // cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2036 // 2037 // // test whether old center is on the band's plane 2038 // if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 2039 // cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 2040 // OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 2041 // } 2042 // radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2043 // if (fabs(radius - CircleRadius) < HULLEPSILON) { 2044 // 2045 // // construct SearchDirection 2046 // SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal); 2047 // helper.CopyVector(&BaseLine->endpoints[0]->node->x); 2048 // for(int i=0;i<3;i++) // just take next different endpoint 2049 // if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) { 2050 // helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x); 2051 // } 2052 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards! 2053 // SearchDirection.Scale(-1.); 2054 // SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2055 // SearchDirection.Normalize(); 2056 // cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2057 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2058 // cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2059 // } 2060 // 2061 // if (LC->SetIndexToVector(&CircleCenter)) { // get cell for the starting atom 2062 // for(int i=0;i<NDIM;i++) // store indices of this cell 2063 // N[i] = LC->n[i]; 2064 // cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2065 // } else { 2066 // cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; 2067 // return; 2068 // } 2069 // // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2070 // cout << Verbose(2) << "LC Intervals:"; 2071 // for (int i=0;i<NDIM;i++) { 2072 // Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2073 // Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2074 // cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2075 // } 2076 // cout << endl; 2077 // for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2078 // for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 2079 // for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2080 // List = LC->GetCurrentCell(); 2081 // cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2082 // if (List != NULL) { 2083 // for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2084 // Candidate = (*Runner); 2085 // 2086 // // check for three unique points 2087 // if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) { 2088 // cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl; 2089 // 2090 // // construct both new centers 2091 // GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)); 2092 // OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2093 // 2094 // if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) { 2095 // helper.CopyVector(&NewNormalVector); 2096 // cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2097 // radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter); 2098 // if (radius < RADIUS*RADIUS) { 2099 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2100 // cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl; 2101 // NewSphereCenter.AddVector(&helper); 2102 // NewSphereCenter.SubtractVector(&CircleCenter); 2103 // cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2104 // 2105 // helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction 2106 // OtherNewSphereCenter.AddVector(&helper); 2107 // OtherNewSphereCenter.SubtractVector(&CircleCenter); 2108 // cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2109 // 2110 // // check both possible centers 2111 // alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection); 2112 // Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection); 2113 // alpha = min(alpha, Otheralpha); 2114 // if (*ShortestAngle > alpha) { 2115 // OptCandidate = Candidate; 2116 // *ShortestAngle = alpha; 2117 // if (alpha != Otheralpha) 2118 // OptCandidateCenter->CopyVector(&NewSphereCenter); 2119 // else 2120 // OptCandidateCenter->CopyVector(&OtherNewSphereCenter); 2121 // cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl; 2122 // } else { 2123 // if (OptCandidate != NULL) 2124 // cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2125 // else 2126 // cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2127 // } 2128 // 2129 // } else { 2130 // cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl; 2131 // } 2132 // } else { 2133 // cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2134 // } 2135 // } else { 2136 // cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl; 2137 // } 2138 // } 2139 // } 2140 // } 2141 // } else { 2142 // cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 2143 // } 2144 // } else { 2145 // cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl; 2146 // } 2147 // 2148 // cout << Verbose(1) << "End of Find_next_suitable_point" << endl; 2149 // }; 2150 1798 /** Checks whether the triangle consisting of the three atoms is already present. 1799 * Searches for the points in Tesselation::PointsOnBoundary and checks their 1800 * lines. If any of the three edges already has two triangles attached, false is 1801 * returned. 1802 * \param *out output stream for debugging 1803 * \param *Candidates endpoints of the triangle candidate 1804 * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two 1805 * triangles exist which is the maximum for three points 1806 */ 1807 int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) { 1808 LineMap::iterator FindLine; 1809 PointMap::iterator FindPoint; 1810 TriangleMap::iterator FindTriangle; 1811 int adjacentTriangleCount = 0; 1812 class BoundaryPointSet *Points[3]; 1813 1814 //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl; 1815 // builds a triangle point set (Points) of the end points 1816 for (int i = 0; i < 3; i++) { 1817 FindPoint = PointsOnBoundary.find(Candidates[i]->nr); 1818 if (FindPoint != PointsOnBoundary.end()) { 1819 Points[i] = FindPoint->second; 1820 } else { 1821 Points[i] = NULL; 1822 } 1823 } 1824 1825 // checks lines between the points in the Points for their adjacent triangles 1826 for (int i = 0; i < 3; i++) { 1827 if (Points[i] != NULL) { 1828 for (int j = i; j < 3; j++) { 1829 if (Points[j] != NULL) { 1830 FindLine = Points[i]->lines.find(Points[j]->node->nr); 1831 if (FindLine != Points[i]->lines.end()) { 1832 for (; FindLine->first == Points[j]->node->nr; FindLine++) { 1833 FindTriangle = FindLine->second->triangles.begin(); 1834 for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 1835 if (( 1836 (FindTriangle->second->endpoints[0] == Points[0]) 1837 || (FindTriangle->second->endpoints[0] == Points[1]) 1838 || (FindTriangle->second->endpoints[0] == Points[2]) 1839 ) && ( 1840 (FindTriangle->second->endpoints[1] == Points[0]) 1841 || (FindTriangle->second->endpoints[1] == Points[1]) 1842 || (FindTriangle->second->endpoints[1] == Points[2]) 1843 ) && ( 1844 (FindTriangle->second->endpoints[2] == Points[0]) 1845 || (FindTriangle->second->endpoints[2] == Points[1]) 1846 || (FindTriangle->second->endpoints[2] == Points[2]) 1847 ) 1848 ) { 1849 adjacentTriangleCount++; 1850 } 1851 } 1852 } 1853 // Only one of the triangle lines must be considered for the triangle count. 1854 *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1855 return adjacentTriangleCount; 1856 1857 } 1858 } 1859 } 1860 } 1861 } 1862 1863 *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1864 return adjacentTriangleCount; 1865 }; 2151 1866 2152 1867 /** This recursive function finds a third point, to form a triangle with two given ones. … … 2177 1892 * @param BaseLine BoundaryLineSet with the current base line 2178 1893 * @param ThirdNode third atom to avoid in search 2179 * @param OptCandidate candidate reference on return 2180 * @param OptCandidateCenter candidate's sphere center on return 1894 * @param candidates list of equally good candidates to return 2181 1895 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate 2182 1896 * @param RADIUS radius of sphere 2183 1897 * @param *LC LinkedCell structure with neighbouring atoms 2184 1898 */ 2185 void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC) 2186 { 1899 void Find_third_point_for_Tesselation( 1900 Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, 1901 class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates, 1902 double *ShortestAngle, const double RADIUS, LinkedCell *LC 1903 ) { 2187 1904 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2188 1905 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 1906 Vector SphereCenter; 2189 1907 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 2190 1908 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 2191 1909 Vector NewNormalVector; // normal vector of the Candidate's triangle 2192 Vector helper ;1910 Vector helper, OptCandidateCenter, OtherOptCandidateCenter; 2193 1911 LinkedAtoms *List = NULL; 2194 1912 double CircleRadius; // radius of this circle … … 2197 1915 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2198 1916 atom *Candidate = NULL; 1917 CandidateForTesselation *optCandidate = NULL; 2199 1918 2200 1919 cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl; 2201 1920 2202 cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;1921 //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2203 1922 2204 1923 // construct center of circle … … 2211 1930 CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x); 2212 1931 2213 // calculate squared radius of circle1932 // calculate squared radius atom *ThirdNode,f circle 2214 1933 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2215 1934 if (radius/4. < RADIUS*RADIUS) { 2216 1935 CircleRadius = RADIUS*RADIUS - radius/4.; 2217 1936 CirclePlaneNormal.Normalize(); 2218 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;1937 //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2219 1938 2220 1939 // test whether old center is on the band's plane … … 2227 1946 2228 1947 // check SearchDirection 2229 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;1948 //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2230 1949 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2231 1950 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; 2232 1951 } 1952 2233 1953 // get cell for the starting atom 2234 1954 if (LC->SetIndexToVector(&CircleCenter)) { 2235 2236 2237 cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;1955 for(int i=0;i<NDIM;i++) // store indices of this cell 1956 N[i] = LC->n[i]; 1957 //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2238 1958 } else { 2239 1959 cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; … … 2241 1961 } 2242 1962 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2243 cout << Verbose(2) << "LC Intervals:";1963 //cout << Verbose(2) << "LC Intervals:"; 2244 1964 for (int i=0;i<NDIM;i++) { 2245 1965 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2246 1966 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2247 cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";1967 //cout << " [" << Nlower[i] << "," << Nupper[i] << "] "; 2248 1968 } 2249 cout << endl;1969 //cout << endl; 2250 1970 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2251 1971 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) … … 2258 1978 2259 1979 // check for three unique points 2260 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {2261 cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;1980 //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl; 1981 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){ 2262 1982 2263 1983 // construct both new centers … … 2265 1985 OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2266 1986 2267 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) { 1987 if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) 1988 && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON) 1989 ) { 2268 1990 helper.CopyVector(&NewNormalVector); 2269 cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;1991 //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2270 1992 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter); 2271 1993 if (radius < RADIUS*RADIUS) { 2272 1994 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2273 cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm()<< "." << endl;1995 //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl; 2274 1996 NewSphereCenter.AddVector(&helper); 2275 1997 NewSphereCenter.SubtractVector(&CircleCenter); 2276 cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2277 2278 helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction 1998 //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 1999 2000 // OtherNewSphereCenter is created by the same vector just in the other direction 2001 helper.Scale(-1.); 2279 2002 OtherNewSphereCenter.AddVector(&helper); 2280 2003 OtherNewSphereCenter.SubtractVector(&CircleCenter); 2281 cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2282 2283 // check both possible centers 2004 //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2005 2284 2006 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2285 2007 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2286 2008 alpha = min(alpha, Otheralpha); 2287 if (*ShortestAngle > alpha) { 2288 OptCandidate = Candidate; 2289 *ShortestAngle = alpha; 2290 if (alpha != Otheralpha) 2291 OptCandidateCenter->CopyVector(&NewSphereCenter); 2292 else 2293 OptCandidateCenter->CopyVector(&OtherNewSphereCenter); 2294 cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl; 2009 // if there is a better candidate, drop the current list and add the new candidate 2010 // otherwise ignore the new candidate and keep the list 2011 if (*ShortestAngle > (alpha - HULLEPSILON)) { 2012 optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter); 2013 if (fabs(alpha - Otheralpha) > MYEPSILON) { 2014 optCandidate->OptCenter.CopyVector(&NewSphereCenter); 2015 optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter); 2016 } else { 2017 optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter); 2018 optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter); 2019 } 2020 // if there is an equal candidate, add it to the list without clearing the list 2021 if ((*ShortestAngle - HULLEPSILON) < alpha) { 2022 candidates->push_back(optCandidate); 2023 cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with " 2024 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2025 } else { 2026 // remove all candidates from the list and then the list itself 2027 class CandidateForTesselation *remover = NULL; 2028 for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) { 2029 remover = *it; 2030 delete(remover); 2031 } 2032 candidates->clear(); 2033 candidates->push_back(optCandidate); 2034 cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with " 2035 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2036 } 2037 *ShortestAngle = alpha; 2038 //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl; 2295 2039 } else { 2296 if (OptCandidate != NULL) 2297 cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2298 else 2299 cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2040 if ((optCandidate != NULL) && (optCandidate->point != NULL)) { 2041 //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl; 2042 } else { 2043 //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2044 } 2300 2045 } 2301 2046 2302 2047 } else { 2303 cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;2048 //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl; 2304 2049 } 2305 2050 } else { 2306 cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;2051 //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2307 2052 } 2308 2053 } else { 2309 if (ThirdNode != NULL) 2310 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2311 else 2312 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl; 2054 if (ThirdNode != NULL) { 2055 //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2056 } else { 2057 //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl; 2058 } 2313 2059 } 2314 2060 } … … 2316 2062 } 2317 2063 } else { 2318 cerr << Verbose( 1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;2064 cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 2319 2065 } 2320 2066 } else { 2321 2067 if (ThirdNode != NULL) 2322 cout << Verbose( 1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;2068 cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 2323 2069 else 2324 cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2325 } 2326 2070 cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2071 } 2072 2073 //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl; 2074 if (candidates->size() > 1) { 2075 candidates->unique(); 2076 candidates->sort(sortCandidates); 2077 } 2327 2078 2328 2079 cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl; 2329 2080 }; 2330 2081 2331 /** Checks whether the triangle consisting of the three atoms is already present. 2332 * Searches for the points in Tesselation::PointsOnBoundary and checks their 2333 * lines. If any of the three edges already has two triangles attached, false is 2334 * returned. 2335 * \param *out output stream for debugging 2336 * \param *a first endpoint 2337 * \param *b second endpoint 2338 * \param *c third endpoint 2339 * \return false - triangle invalid due to edge criteria, true - triangle may be added. 2340 */ 2341 bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c) { 2342 LineMap::iterator TryAndError; 2343 PointTestPair InsertPoint; 2344 bool Present[3]; 2345 2346 *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl; 2347 TPS[0] = new class BoundaryPointSet(a); 2348 TPS[1] = new class BoundaryPointSet(b); 2349 TPS[2] = new class BoundaryPointSet(c); 2350 for (int i=0;i<3;i++) { // check through all endpoints 2351 InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i])); 2352 Present[i] = !InsertPoint.second; 2353 if (Present[i]) { // if new point was not present before, increase counter 2354 delete TPS[i]; 2355 *out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl; 2356 TPS[i] = (InsertPoint.first)->second; 2357 } 2358 } 2359 2360 // check lines 2361 for (int i=0;i<3;i++) 2362 if (Present[i]) 2363 for (int j=i;j<3;j++) 2364 if (Present[j]) { 2365 TryAndError = TPS[i]->lines.find(TPS[j]->node->nr); 2366 if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) { 2367 *out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl; 2368 *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 2369 return false; 2370 } 2082 struct Intersection { 2083 Vector x1; 2084 Vector x2; 2085 Vector x3; 2086 Vector x4; 2087 }; 2088 2089 /** 2090 * Intersection calculation function. 2091 * 2092 * @param x to find the result for 2093 * @param function parameter 2094 */ 2095 double MinIntersectDistance(const gsl_vector * x, void *params) { 2096 double retval = 0; 2097 struct Intersection *I = (struct Intersection *)params; 2098 Vector intersection; 2099 Vector SideA,SideB,HeightA, HeightB; 2100 for (int i=0;i<NDIM;i++) 2101 intersection.x[i] = gsl_vector_get(x, i); 2102 2103 SideA.CopyVector(&(I->x1)); 2104 SideA.SubtractVector(&I->x2); 2105 HeightA.CopyVector(&intersection); 2106 HeightA.SubtractVector(&I->x1); 2107 HeightA.ProjectOntoPlane(&SideA); 2108 2109 SideB.CopyVector(&I->x3); 2110 SideB.SubtractVector(&I->x4); 2111 HeightB.CopyVector(&intersection); 2112 HeightB.SubtractVector(&I->x3); 2113 HeightB.ProjectOntoPlane(&SideB); 2114 2115 retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB); 2116 //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl; 2117 2118 return retval; 2119 }; 2120 2121 2122 /** 2123 * Calculates whether there is an intersection between two lines. The first line 2124 * always goes through point 1 and point 2 and the second line is given by the 2125 * connection between point 4 and point 5. 2126 * 2127 * @param point 1 of line 1 2128 * @param point 2 of line 1 2129 * @param point 1 of line 2 2130 * @param point 2 of line 2 2131 * 2132 * @return true if there is an intersection between the given lines, false otherwise 2133 */ 2134 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) { 2135 bool result; 2136 2137 struct Intersection par; 2138 par.x1.CopyVector(&point1); 2139 par.x2.CopyVector(&point2); 2140 par.x3.CopyVector(&point3); 2141 par.x4.CopyVector(&point4); 2142 2143 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; 2144 gsl_multimin_fminimizer *s = NULL; 2145 gsl_vector *ss, *x; 2146 gsl_multimin_function minex_func; 2147 2148 size_t iter = 0; 2149 int status; 2150 double size; 2151 2152 /* Starting point */ 2153 x = gsl_vector_alloc(NDIM); 2154 gsl_vector_set(x, 0, point1.x[0]); 2155 gsl_vector_set(x, 1, point1.x[1]); 2156 gsl_vector_set(x, 2, point1.x[2]); 2157 2158 /* Set initial step sizes to 1 */ 2159 ss = gsl_vector_alloc(NDIM); 2160 gsl_vector_set_all(ss, 1.0); 2161 2162 /* Initialize method and iterate */ 2163 minex_func.n = NDIM; 2164 minex_func.f = &MinIntersectDistance; 2165 minex_func.params = (void *)∥ 2166 2167 s = gsl_multimin_fminimizer_alloc(T, NDIM); 2168 gsl_multimin_fminimizer_set(s, &minex_func, x, ss); 2169 2170 do { 2171 iter++; 2172 status = gsl_multimin_fminimizer_iterate(s); 2173 2174 if (status) { 2175 break; 2371 2176 } 2372 *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 2373 return true; 2374 }; 2375 2177 2178 size = gsl_multimin_fminimizer_size(s); 2179 status = gsl_multimin_test_size(size, 1e-2); 2180 2181 if (status == GSL_SUCCESS) { 2182 cout << Verbose(2) << "converged to minimum" << endl; 2183 } 2184 } while (status == GSL_CONTINUE && iter < 100); 2185 2186 // check whether intersection is in between or not 2187 Vector intersection, SideA, SideB, HeightA, HeightB; 2188 double t1, t2; 2189 for (int i = 0; i < NDIM; i++) { 2190 intersection.x[i] = gsl_vector_get(s->x, i); 2191 } 2192 2193 SideA.CopyVector(&par.x2); 2194 SideA.SubtractVector(&par.x1); 2195 HeightA.CopyVector(&intersection); 2196 HeightA.SubtractVector(&par.x1); 2197 2198 t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA); 2199 2200 SideB.CopyVector(&par.x4); 2201 SideB.SubtractVector(&par.x3); 2202 HeightB.CopyVector(&intersection); 2203 HeightB.SubtractVector(&par.x3); 2204 2205 t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB); 2206 2207 cout << Verbose(2) << "Intersection " << intersection << " is at " 2208 << t1 << " for (" << point1 << "," << point2 << ") and at " 2209 << t2 << " for (" << point3 << "," << point4 << "): "; 2210 2211 if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) { 2212 cout << "true intersection." << endl; 2213 result = true; 2214 } else { 2215 cout << "intersection out of region of interest." << endl; 2216 result = false; 2217 } 2218 2219 // free minimizer stuff 2220 gsl_vector_free(x); 2221 gsl_vector_free(ss); 2222 gsl_multimin_fminimizer_free(s); 2223 2224 return result; 2225 } 2376 2226 2377 2227 /** Finds the second point of starting triangle. … … 2400 2250 } 2401 2251 // then go through the current and all neighbouring cells and check the contained atoms for possible candidates 2402 cout << Verbose(2) << "LC Intervals:"; 2252 cout << Verbose(3) << "LC Intervals from ["; 2253 for (int i=0;i<NDIM;i++) { 2254 cout << " " << N[i] << "<->" << LC->N[i]; 2255 } 2256 cout << "] :"; 2403 2257 for (int i=0;i<NDIM;i++) { 2404 2258 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; … … 2407 2261 } 2408 2262 cout << endl; 2263 2264 2409 2265 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 2410 2266 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) … … 2415 2271 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 2416 2272 Candidate = (*Runner); 2417 // check if we only have one unique point yet ... 2418 if (a != Candidate) { 2419 cout << Verbose(3) << "Current candidate is " << *Candidate << ": "; 2420 AngleCheck.CopyVector(&(Candidate->x)); 2421 AngleCheck.SubtractVector(&(a->x)); 2422 norm = AngleCheck.Norm(); 2273 // check if we only have one unique point yet ... 2274 if (a != Candidate) { 2275 // Calculate center of the circle with radius RADIUS through points a and Candidate 2276 Vector OrthogonalizedOben, a_Candidate, Center; 2277 double distance, scaleFactor; 2278 2279 OrthogonalizedOben.CopyVector(&Oben); 2280 a_Candidate.CopyVector(&(a->x)); 2281 a_Candidate.SubtractVector(&(Candidate->x)); 2282 OrthogonalizedOben.ProjectOntoPlane(&a_Candidate); 2283 OrthogonalizedOben.Normalize(); 2284 distance = 0.5 * a_Candidate.Norm(); 2285 scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance))); 2286 OrthogonalizedOben.Scale(scaleFactor); 2287 2288 Center.CopyVector(&(Candidate->x)); 2289 Center.AddVector(&(a->x)); 2290 Center.Scale(0.5); 2291 Center.AddVector(&OrthogonalizedOben); 2292 2293 AngleCheck.CopyVector(&Center); 2294 AngleCheck.SubtractVector(&(a->x)); 2295 norm = a_Candidate.Norm(); 2423 2296 // second point shall have smallest angle with respect to Oben vector 2424 if (norm < RADIUS ) {2297 if (norm < RADIUS*2.) { 2425 2298 angle = AngleCheck.Angle(&Oben); 2426 2299 if (angle < Storage[0]) { 2427 //cout << Verbose( 1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);2428 cout << "Is a better candidate with distance " << norm << " and " << angle<< ".\n";2300 //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 2301 cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n"; 2429 2302 Opt_Candidate = Candidate; 2430 Storage[0] = AngleCheck.Angle(&Oben);2431 //cout << Verbose( 1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);2303 Storage[0] = angle; 2304 //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]); 2432 2305 } else { 2433 cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;2306 //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl; 2434 2307 } 2435 2308 } else { 2436 cout << "Refused due to Radius " << norm << endl;2309 //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl; 2437 2310 } 2311 } else { 2312 //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl; 2438 2313 } 2439 2314 } 2315 } else { 2316 cout << Verbose(3) << "Linked cell list is empty." << endl; 2440 2317 } 2441 2318 } … … 2455 2332 int i = 0; 2456 2333 LinkedAtoms *List = NULL; 2457 atom* FirstPoint ;2458 atom* SecondPoint ;2334 atom* FirstPoint = NULL; 2335 atom* SecondPoint = NULL; 2459 2336 atom* MaxAtom[NDIM]; 2460 2337 double max_coordinate[NDIM]; … … 2463 2340 Vector Chord; 2464 2341 Vector SearchDirection; 2465 Vector OptCandidateCenter;2466 2342 2467 2343 Oben.Zero(); … … 2481 2357 if (List != NULL) { 2482 2358 for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) { 2483 cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;2484 2359 if ((*Runner)->x.x[i] > max_coordinate[i]) { 2360 cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl; 2485 2361 max_coordinate[i] = (*Runner)->x.x[i]; 2486 2362 MaxAtom[i] = (*Runner); … … 2497 2373 cout << i << ": " << *MaxAtom[i] << "\t"; 2498 2374 cout << endl; 2499 const int k = 1; // arbitrary choice 2500 Oben.x[k] = 1.; 2501 FirstPoint = MaxAtom[k]; 2502 cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl; 2503 2504 // add first point 2505 AddTrianglePoint(FirstPoint, 0); 2506 2507 double ShortestAngle; 2508 atom* Opt_Candidate = NULL; 2509 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant. 2510 2511 Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 2512 SecondPoint = Opt_Candidate; 2513 cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n"; 2514 2515 // add second point and first baseline 2516 AddTrianglePoint(SecondPoint, 1); 2517 AddTriangleLine(TPS[0], TPS[1], 0); 2518 2519 helper.CopyVector(&(FirstPoint->x)); 2520 helper.SubtractVector(&(SecondPoint->x)); 2521 helper.Normalize(); 2522 Oben.ProjectOntoPlane(&helper); 2523 Oben.Normalize(); 2524 helper.VectorProduct(&Oben); 2525 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2526 2527 Chord.CopyVector(&(FirstPoint->x)); // bring into calling function 2528 Chord.SubtractVector(&(SecondPoint->x)); 2529 double radius = Chord.ScalarProduct(&Chord); 2530 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 2531 helper.CopyVector(&Oben); 2532 helper.Scale(CircleRadius); 2533 // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized) 2534 2535 cout << Verbose(2) << "Looking for third point candidates \n"; 2536 // look in one direction of baseline for initial candidate 2537 Opt_Candidate = NULL; 2538 SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ... 2539 2540 cout << Verbose(1) << "Looking for third point candidates ...\n"; 2541 Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC); 2542 cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl; 2543 2544 // add third point 2545 AddTrianglePoint(Opt_Candidate, 2); 2546 2547 // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate 2548 2549 // Finally, we only have to add the found further lines 2550 AddTriangleLine(TPS[1], TPS[2], 1); 2551 AddTriangleLine(TPS[0], TPS[2], 2); 2552 // ... and triangles to the Maps of the Tesselation class 2553 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2554 AddTriangleToLines(); 2555 // ... and calculate its normal vector (with correct orientation) 2556 OptCandidateCenter.Scale(-1.); 2557 cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl; 2558 BTS->GetNormalVector(OptCandidateCenter); 2559 cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n"; 2560 cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl; 2375 2376 BTS = NULL; 2377 CandidateList *Opt_Candidates = new CandidateList(); 2378 for (int k=0;k<NDIM;k++) { 2379 Oben.x[k] = 1.; 2380 FirstPoint = MaxAtom[k]; 2381 cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl; 2382 2383 double ShortestAngle; 2384 atom* Opt_Candidate = NULL; 2385 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant. 2386 2387 Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 2388 SecondPoint = Opt_Candidate; 2389 if (SecondPoint == NULL) // have we found a second point? 2390 continue; 2391 else 2392 cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n"; 2393 2394 helper.CopyVector(&(FirstPoint->x)); 2395 helper.SubtractVector(&(SecondPoint->x)); 2396 helper.Normalize(); 2397 Oben.ProjectOntoPlane(&helper); 2398 Oben.Normalize(); 2399 helper.VectorProduct(&Oben); 2400 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2401 2402 Chord.CopyVector(&(FirstPoint->x)); // bring into calling function 2403 Chord.SubtractVector(&(SecondPoint->x)); 2404 double radius = Chord.ScalarProduct(&Chord); 2405 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 2406 helper.CopyVector(&Oben); 2407 helper.Scale(CircleRadius); 2408 // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized) 2409 2410 // look in one direction of baseline for initial candidate 2411 SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ... 2412 2413 // adding point 1 and point 2 and the line between them 2414 AddTrianglePoint(FirstPoint, 0); 2415 AddTrianglePoint(SecondPoint, 1); 2416 AddTriangleLine(TPS[0], TPS[1], 0); 2417 2418 //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n"; 2419 Find_third_point_for_Tesselation( 2420 Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC 2421 ); 2422 cout << Verbose(1) << "List of third Points is "; 2423 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2424 cout << " " << *(*it)->point; 2425 } 2426 cout << endl; 2427 2428 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2429 // add third triangle point 2430 AddTrianglePoint((*it)->point, 2); 2431 // add the second and third line 2432 AddTriangleLine(TPS[1], TPS[2], 1); 2433 AddTriangleLine(TPS[0], TPS[2], 2); 2434 // ... and triangles to the Maps of the Tesselation class 2435 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2436 AddTriangle(); 2437 // ... and calculate its normal vector (with correct orientation) 2438 (*it)->OptCenter.Scale(-1.); 2439 cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl; 2440 BTS->GetNormalVector((*it)->OptCenter); // vector to compare with should point inwards 2441 cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " 2442 << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n"; 2443 2444 // if we do not reach the end with the next step of iteration, we need to setup a new first line 2445 if (it != Opt_Candidates->end()--) { 2446 FirstPoint = (*it)->BaseLine->endpoints[0]->node; 2447 SecondPoint = (*it)->point; 2448 // adding point 1 and point 2 and the line between them 2449 AddTrianglePoint(FirstPoint, 0); 2450 AddTrianglePoint(SecondPoint, 1); 2451 AddTriangleLine(TPS[0], TPS[1], 0); 2452 } 2453 cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl; 2454 } 2455 if (BTS != NULL) // we have created one starting triangle 2456 break; 2457 else { 2458 // remove all candidates from the list and then the list itself 2459 class CandidateForTesselation *remover = NULL; 2460 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2461 remover = *it; 2462 delete(remover); 2463 } 2464 Opt_Candidates->clear(); 2465 } 2466 } 2467 2468 // remove all candidates from the list and then the list itself 2469 class CandidateForTesselation *remover = NULL; 2470 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2471 remover = *it; 2472 delete(remover); 2473 } 2474 delete(Opt_Candidates); 2561 2475 cout << Verbose(1) << "End of Find_starting_triangle\n"; 2562 2476 }; 2477 2478 /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected. 2479 * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not 2480 * make it bigger (i.e. closing one (the baseline) and opening two new ones). 2481 * \param TPS[3] nodes of the triangle 2482 * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create) 2483 */ 2484 bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3]) 2485 { 2486 bool result = false; 2487 int counter = 0; 2488 2489 // check all three points 2490 for (int i=0;i<3;i++) 2491 for (int j=i+1; j<3; j++) { 2492 if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) { // there already is a line 2493 LineMap::iterator FindLine; 2494 pair<LineMap::iterator,LineMap::iterator> FindPair; 2495 FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr); 2496 for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) { 2497 // If there is a line with less than two attached triangles, we don't need a new line. 2498 if (FindLine->second->TrianglesCount < 2) { 2499 counter++; 2500 break; // increase counter only once per edge 2501 } 2502 } 2503 } else { // no line 2504 cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl; 2505 result = true; 2506 } 2507 } 2508 if (counter > 1) { 2509 cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl; 2510 result = true; 2511 } 2512 return result; 2513 }; 2514 2563 2515 2564 2516 /** This function finds a triangle to a line, adjacent to an existing one. … … 2576 2528 const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC) 2577 2529 { 2578 cout << Verbose( 1) << "Begin of Find_next_suitable_triangle\n";2530 cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n"; 2579 2531 ofstream *tempstream = NULL; 2580 2532 char NumberName[255]; 2581 2582 atom* Opt_Candidate = NULL; 2583 Vector OptCandidateCenter; 2533 bool result = true; 2534 CandidateList *Opt_Candidates = new CandidateList(); 2584 2535 2585 2536 Vector CircleCenter; … … 2589 2540 Vector helper; 2590 2541 atom *ThirdNode = NULL; 2542 LineMap::iterator testline; 2591 2543 double ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2592 2544 double radius, CircleRadius; … … 2620 2572 OldSphereCenter.AddVector(&helper); 2621 2573 OldSphereCenter.SubtractVector(&CircleCenter); 2622 cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;2574 //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2623 2575 2624 2576 // construct SearchDirection … … 2626 2578 helper.CopyVector(&Line.endpoints[0]->node->x); 2627 2579 helper.SubtractVector(&ThirdNode->x); 2628 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) 2580 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 2629 2581 SearchDirection.Scale(-1.); 2630 2582 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2631 2583 SearchDirection.Normalize(); 2632 2584 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2633 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2585 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2586 // rotated the wrong way! 2634 2587 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2635 2588 } 2636 2589 2637 2590 // add third point 2638 cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl; 2639 Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC); 2591 Find_third_point_for_Tesselation( 2592 T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates, 2593 &ShortestAngle, RADIUS, LC 2594 ); 2640 2595 2641 2596 } else { … … 2643 2598 } 2644 2599 2645 if (Opt_Candidate == NULL) {2600 if (Opt_Candidates->begin() == Opt_Candidates->end()) { 2646 2601 cerr << "WARNING: Could not find a suitable candidate." << endl; 2647 2602 return false; 2648 2603 } 2649 cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl; 2650 2651 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2652 bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node); 2653 2654 if (flag) { // if so, add 2655 AddTrianglePoint(Opt_Candidate, 0); 2656 AddTrianglePoint(Line.endpoints[0]->node, 1); 2657 AddTrianglePoint(Line.endpoints[1]->node, 2); 2658 2659 AddTriangleLine(TPS[0], TPS[1], 0); 2660 AddTriangleLine(TPS[0], TPS[2], 1); 2661 AddTriangleLine(TPS[1], TPS[2], 2); 2662 2663 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2664 AddTriangleToLines(); 2665 2666 OptCandidateCenter.Scale(-1.); 2667 BTS->GetNormalVector(OptCandidateCenter); 2668 OptCandidateCenter.Scale(-1.); 2669 2670 cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl; 2671 cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl; 2672 } else { // else, yell and do nothing 2673 cout << Verbose(1) << "This triangle consisting of "; 2674 cout << *Opt_Candidate << ", "; 2675 cout << *Line.endpoints[0]->node << " and "; 2676 cout << *Line.endpoints[1]->node << " "; 2677 cout << "is invalid!" << endl; 2678 return false; 2679 } 2680 2681 if (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration 2682 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name); 2683 if (DoTecplotOutput) { 2684 string NameofTempFile(tempbasename); 2685 NameofTempFile.append(NumberName); 2686 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos)) 2604 cout << Verbose(1) << "Third Points are "; 2605 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2606 cout << " " << *(*it)->point; 2607 } 2608 cout << endl; 2609 2610 BoundaryLineSet *BaseRay = &Line; 2611 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2612 cout << Verbose(1) << " Third point candidate is " << *(*it)->point 2613 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2614 cout << Verbose(1) << " Baseline is " << *BaseRay << endl; 2615 2616 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2617 atom *AtomCandidates[3]; 2618 AtomCandidates[0] = (*it)->point; 2619 AtomCandidates[1] = BaseRay->endpoints[0]->node; 2620 AtomCandidates[2] = BaseRay->endpoints[1]->node; 2621 int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates); 2622 2623 BTS = NULL; 2624 // If there is no triangle, add it regularly. 2625 if (existentTrianglesCount == 0) { 2626 AddTrianglePoint((*it)->point, 0); 2627 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); 2628 AddTrianglePoint(BaseRay->endpoints[1]->node, 2); 2629 2630 AddTriangleLine(TPS[0], TPS[1], 0); 2631 AddTriangleLine(TPS[0], TPS[2], 1); 2632 AddTriangleLine(TPS[1], TPS[2], 2); 2633 2634 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2635 AddTriangle(); 2636 (*it)->OptCenter.Scale(-1.); 2637 BTS->GetNormalVector((*it)->OptCenter); 2638 (*it)->OptCenter.Scale(-1.); 2639 2640 cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector 2641 << " for this triangle ... " << endl; 2642 //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl; 2643 } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time. 2644 AddTrianglePoint((*it)->point, 0); 2645 AddTrianglePoint(BaseRay->endpoints[0]->node, 1); 2646 AddTrianglePoint(BaseRay->endpoints[1]->node, 2); 2647 2648 // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1) 2649 // i.e. at least one of the three lines must be present with TriangleCount <= 1 2650 if (CheckLineCriteriaforDegeneratedTriangle(TPS)) { 2651 AddTriangleLine(TPS[0], TPS[1], 0); 2652 AddTriangleLine(TPS[0], TPS[2], 1); 2653 AddTriangleLine(TPS[1], TPS[2], 2); 2654 2655 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2656 AddTriangle(); 2657 2658 (*it)->OtherOptCenter.Scale(-1.); 2659 BTS->GetNormalVector((*it)->OtherOptCenter); 2660 (*it)->OtherOptCenter.Scale(-1.); 2661 2662 cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector 2663 << " for this triangle ... " << endl; 2664 cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl; 2665 } else { 2666 cout << Verbose(1) << "WARNING: This triangle consisting of "; 2667 cout << *(*it)->point << ", "; 2668 cout << *BaseRay->endpoints[0]->node << " and "; 2669 cout << *BaseRay->endpoints[1]->node << " "; 2670 cout << "exists and is not added, as it does not seem helpful!" << endl; 2671 result = false; 2672 } 2673 } else { 2674 cout << Verbose(1) << "This triangle consisting of "; 2675 cout << *(*it)->point << ", "; 2676 cout << *BaseRay->endpoints[0]->node << " and "; 2677 cout << *BaseRay->endpoints[1]->node << " "; 2678 cout << "is invalid!" << endl; 2679 result = false; 2680 } 2681 2682 if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration 2683 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name); 2684 if (DoTecplotOutput) { 2685 string NameofTempFile(tempbasename); 2686 NameofTempFile.append(NumberName); 2687 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos)) 2687 2688 NameofTempFile.erase(npos, 1); 2688 NameofTempFile.append(TecplotSuffix);2689 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";2690 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);2691 write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);2692 tempstream->close();2693 tempstream->flush();2694 delete(tempstream);2695 }2696 2697 if (DoRaster3DOutput) {2698 string NameofTempFile(tempbasename);2699 NameofTempFile.append(NumberName);2700 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))2689 NameofTempFile.append(TecplotSuffix); 2690 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 2691 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 2692 write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten); 2693 tempstream->close(); 2694 tempstream->flush(); 2695 delete(tempstream); 2696 } 2697 2698 if (DoRaster3DOutput) { 2699 string NameofTempFile(tempbasename); 2700 NameofTempFile.append(NumberName); 2701 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos)) 2701 2702 NameofTempFile.erase(npos, 1); 2702 NameofTempFile.append(Raster3DSuffix); 2703 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 2704 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 2705 write_raster3d_file(out, tempstream, this, mol); 2706 // include the current position of the virtual sphere in the temporary raster3d file 2707 // make the circumsphere's center absolute again 2708 helper.CopyVector(&Line.endpoints[0]->node->x); 2709 helper.AddVector(&Line.endpoints[1]->node->x); 2710 helper.Scale(0.5); 2711 OptCandidateCenter.AddVector(&helper); 2712 Vector *center = mol->DetermineCenterOfAll(out); 2713 OptCandidateCenter.AddVector(center); 2714 delete(center); 2715 // and add to file plus translucency object 2716 *tempstream << "# current virtual sphere\n"; 2717 *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 2718 *tempstream << "2\n " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n"; 2719 *tempstream << "9\n terminating special property\n"; 2720 tempstream->close(); 2721 tempstream->flush(); 2722 delete(tempstream); 2723 } 2724 if (DoTecplotOutput || DoRaster3DOutput) 2725 TriangleFilesWritten++; 2726 } 2727 2728 cout << Verbose(1) << "End of Find_next_suitable_triangle\n"; 2729 return true; 2703 NameofTempFile.append(Raster3DSuffix); 2704 cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 2705 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 2706 write_raster3d_file(out, tempstream, this, mol); 2707 // include the current position of the virtual sphere in the temporary raster3d file 2708 // make the circumsphere's center absolute again 2709 helper.CopyVector(&BaseRay->endpoints[0]->node->x); 2710 helper.AddVector(&BaseRay->endpoints[1]->node->x); 2711 helper.Scale(0.5); 2712 (*it)->OptCenter.AddVector(&helper); 2713 Vector *center = mol->DetermineCenterOfAll(out); 2714 (*it)->OptCenter.AddVector(center); 2715 delete(center); 2716 // and add to file plus translucency object 2717 *tempstream << "# current virtual sphere\n"; 2718 *tempstream << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 2719 *tempstream << "2\n " << (*it)->OptCenter.x[0] << " " 2720 << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2] 2721 << "\t" << RADIUS << "\t1 0 0\n"; 2722 *tempstream << "9\n terminating special property\n"; 2723 tempstream->close(); 2724 tempstream->flush(); 2725 delete(tempstream); 2726 } 2727 if (DoTecplotOutput || DoRaster3DOutput) 2728 TriangleFilesWritten++; 2729 } 2730 2731 // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2732 BaseRay = BLS[0]; 2733 } 2734 2735 // remove all candidates from the list and then the list itself 2736 class CandidateForTesselation *remover = NULL; 2737 for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) { 2738 remover = *it; 2739 delete(remover); 2740 } 2741 delete(Opt_Candidates); 2742 cout << Verbose(0) << "End of Find_next_suitable_triangle\n"; 2743 return result; 2730 2744 }; 2745 2746 /** 2747 * Sort function for the candidate list. 2748 */ 2749 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) { 2750 Vector BaseLineVector, OrthogonalVector, helper; 2751 if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 2752 cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 2753 //return false; 2754 exit(1); 2755 } 2756 // create baseline vector 2757 BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x)); 2758 BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2759 BaseLineVector.Normalize(); 2760 2761 // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 2762 helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2763 helper.SubtractVector(&(candidate1->point->x)); 2764 OrthogonalVector.CopyVector(&helper); 2765 helper.VectorProduct(&BaseLineVector); 2766 OrthogonalVector.SubtractVector(&helper); 2767 OrthogonalVector.Normalize(); 2768 2769 // calculate both angles and correct with in-plane vector 2770 helper.CopyVector(&(candidate1->point->x)); 2771 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2772 double phi = BaseLineVector.Angle(&helper); 2773 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 2774 phi = 2.*M_PI - phi; 2775 } 2776 helper.CopyVector(&(candidate2->point->x)); 2777 helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x)); 2778 double psi = BaseLineVector.Angle(&helper); 2779 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 2780 psi = 2.*M_PI - psi; 2781 } 2782 2783 cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl; 2784 cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl; 2785 2786 // return comparison 2787 return phi < psi; 2788 } 2731 2789 2732 2790 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule. … … 2734 2792 * \param *mol molecule structure with Atom's and Bond's 2735 2793 * \param *Tess Tesselation filled with points, lines and triangles on boundary on return 2736 * \param *LCList linked cell list of all atoms2737 2794 * \param *filename filename prefix for output of vertex data 2738 2795 * \para RADIUS radius of the virtual sphere … … 2742 2799 int N = 0; 2743 2800 bool freeTess = false; 2801 bool freeLC = false; 2744 2802 *out << Verbose(1) << "Entering search for non convex hull. " << endl; 2745 2803 if (Tess == NULL) { … … 2748 2806 freeTess = true; 2749 2807 } 2750 bool freeLC = false;2751 2808 LineMap::iterator baseline; 2809 LineMap::iterator testline; 2752 2810 *out << Verbose(0) << "Begin of Find_non_convex_border\n"; 2753 2811 bool flag = false; // marks whether we went once through all baselines without finding any without two triangles … … 2769 2827 cerr << "WARNING: Find_next_suitable_triangle failed." << endl; 2770 2828 } else { 2771 cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl; 2772 } 2829 //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl; 2830 if (baseline->second->TrianglesCount != 2) 2831 cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl; 2832 } 2833 2773 2834 N++; 2774 2835 baseline++; … … 2797 2858 } 2798 2859 } else { 2799 cerr << "ERROR: Could definately not find all necessary triangles!" << endl; 2800 } 2860 cerr << "ERROR: Could definitively not find all necessary triangles!" << endl; 2861 } 2862 2863 cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl; 2864 int counter = 0; 2865 for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) { 2866 if (testline->second->TrianglesCount != 2) { 2867 cout << Verbose(2) << *testline->second << "\t" << testline->second->TrianglesCount << endl; 2868 counter++; 2869 } 2870 } 2871 if (counter == 0) 2872 cout << Verbose(2) << "None." << endl; 2873 2801 2874 if (freeTess) 2802 2875 delete(Tess); -
src/boundary.hpp
r51c910 r36ec71 5 5 class BoundaryLineSet; 6 6 class BoundaryTriangleSet; 7 class CandidateForTesselation; 7 8 8 9 // include config.h … … 17 18 18 19 #include <gsl/gsl_poly.h> 19 #include <gsl/gsl_matrix.h>20 #include <gsl/gsl_linalg.h>21 #include <gsl/gsl_multimin.h>22 #include <gsl/gsl_permutation.h>23 20 24 21 #include "linkedcell.hpp" … … 78 75 }; 79 76 77 78 class CandidateForTesselation { 79 public : 80 CandidateForTesselation(atom* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter); 81 ~CandidateForTesselation(); 82 atom *point; 83 BoundaryLineSet *BaseLine; 84 Vector OptCenter; 85 Vector OtherOptCenter; 86 }; 87 88 80 89 class Tesselation { 81 90 public: … … 89 98 void AddTrianglePoint(atom* Candidate, int n); 90 99 void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n); 91 void AddTriangleToLines(); 100 void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n); 101 void AddTriangle(); 92 102 void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC); 93 103 bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC); 94 bool CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c);104 int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]); 95 105 void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent, int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol); 96 106 … … 117 127 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem); 118 128 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity); 119 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC List, const char *tempbasename, const double RADIUS);129 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS); 120 130 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC); 131 bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess); 132 bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4); 133 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2); 121 134 122 135 #endif /*BOUNDARY_HPP_*/ -
src/builder.cpp
r51c910 r36ec71 307 307 cin >> x.x[i]; 308 308 } 309 // center310 mol->CenterInBox((ofstream *)&cout, &x);311 309 // update Box of atoms by boundary 312 310 mol->SetBoxDimension(&x); 311 // center 312 mol->CenterInBox((ofstream *)&cout); 313 313 break; 314 314 } … … 753 753 static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration) 754 754 { 755 atom *first ;755 atom *first = NULL; 756 756 Vector x,y,z,n; // coordinates for absolute point in cell volume 757 757 int j, axis, count, faktor; … … 1305 1305 * \return exit code (0 - successful, all else - something's wrong) 1306 1306 */ 1307 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char * ConfigFileName, char *&PathToDatabases)1307 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName) 1308 1308 { 1309 1309 Vector x,y,z,n; // coordinates for absolute point in cell volume … … 1320 1320 clock_t start,end; 1321 1321 int argptr; 1322 PathToDatabases = LocalPath;1322 strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1); 1323 1323 1324 1324 // simply create a new molecule, wherein the config file is loaded and the manipulation takes place … … 1343 1343 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl; 1344 1344 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl; 1345 cout << "\t-b x 1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;1345 cout << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl; 1346 1346 cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl; 1347 1347 cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl; … … 1363 1363 cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl; 1364 1364 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; 1365 cout << "\t-T <file> Store temperatures from the config file in <file>." << endl; 1365 cout << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl; 1366 cout << "\t-S <file> Store temperatures from the config file in <file>." << endl; 1366 1367 cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl; 1367 1368 cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; … … 1385 1386 } else { 1386 1387 cout << "Using " << argv[argptr] << " as elements database." << endl; 1387 PathToDatabases = argv[argptr];1388 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1); 1388 1389 argptr+=1; 1389 1390 } … … 1402 1403 1403 1404 // 2. Parse the element database 1404 if (periode->LoadPeriodentafel( PathToDatabases)) {1405 if (periode->LoadPeriodentafel(configuration.databasepath)) { 1405 1406 cout << Verbose(0) << "Element list loaded successfully." << endl; 1406 1407 //periode->Output((ofstream *)&cout); … … 1421 1422 } else { 1422 1423 cout << "Empty configuration file." << endl; 1423 strcpy(ConfigFileName, argv[1]);1424 ConfigFileName = argv[1]; 1424 1425 config_present = empty; 1425 1426 output.close(); … … 1427 1428 } else { 1428 1429 test.close(); 1429 strcpy(ConfigFileName, argv[1]);1430 ConfigFileName = argv[1]; 1430 1431 cout << Verbose(1) << "Specified config file found, parsing ... "; 1431 1432 switch (configuration.TestSyntax(ConfigFileName, periode, mol)) { … … 1576 1577 } else { 1577 1578 class Tesselation T; 1578 int N = 15;1579 int number = 100;1580 1579 string filename(argv[argptr+1]); 1581 1580 filename.append(".csv"); 1582 1581 cout << Verbose(0) << "Evaluating non-convex envelope."; 1583 1582 cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl; 1584 LinkedCell LCList(mol, atof(argv[argptr]) ); // \NOTE not twice the radius??1583 LinkedCell LCList(mol, atof(argv[argptr])*2.); 1585 1584 Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr])); 1586 FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());1585 //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str()); 1587 1586 argptr+=2; 1588 1587 } 1589 1588 break; 1590 case ' T':1589 case 'S': 1591 1590 ExitFlag = 1; 1592 1591 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 1593 1592 ExitFlag = 255; 1594 cerr << "Not enough or invalid arguments given for storing tempature: - T<temperature file>" << endl;1593 cerr << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl; 1595 1594 } else { 1596 1595 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl; … … 1667 1666 x.x[i] = atof(argv[argptr+i]); 1668 1667 mol->Translate((const Vector *)&x); 1668 argptr+=3; 1669 } 1670 case 'T': 1671 ExitFlag = 1; 1672 if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) { 1673 ExitFlag = 255; 1674 cerr << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl; 1675 } else { 1676 ExitFlag = 1; 1677 SaveFlag = true; 1678 cout << Verbose(1) << "Translating all ions periodically to new origin." << endl; 1679 for (int i=NDIM;i--;) 1680 x.x[i] = atof(argv[argptr+i]); 1681 mol->TranslatePeriodically((const Vector *)&x); 1669 1682 argptr+=3; 1670 1683 } … … 1699 1712 case 'b': 1700 1713 ExitFlag = 1; 1701 if ((argptr+ 2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {1714 if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) { 1702 1715 ExitFlag = 255; 1703 cerr << "Not enough or invalid arguments given for centering in box: -b < length_x> <length_y> <length_z>" << endl;1716 cerr << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl; 1704 1717 } else { 1705 1718 SaveFlag = true; 1706 1719 j = -1; 1707 1720 cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl; 1708 j=-1; 1709 for (int i=0;i<NDIM;i++) { 1710 j += i+1; 1711 x.x[i] = atof(argv[argptr++]); 1712 mol->cell_size[j] += x.x[i]*2.; 1721 for (int i=0;i<6;i++) { 1722 mol->cell_size[i] = atof(argv[argptr+i]); 1713 1723 } 1714 1724 // center 1715 mol->CenterInBox((ofstream *)&cout, &x); 1716 // update Box of atoms by boundary 1717 mol->SetBoxDimension(&x); 1725 mol->CenterInBox((ofstream *)&cout); 1726 argptr+=6; 1718 1727 } 1719 1728 break; … … 1735 1744 for (int i=0;i<NDIM;i++) { 1736 1745 j += i+1; 1737 x.x[i] = atof(argv[argptr+ +]);1746 x.x[i] = atof(argv[argptr+i]); 1738 1747 mol->cell_size[j] += x.x[i]*2.; 1739 1748 } 1740 1749 mol->Translate((const Vector *)&x); 1750 argptr+=3; 1741 1751 } 1742 1752 break; … … 1745 1755 SaveFlag = true; 1746 1756 cout << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl; 1757 x.Zero(); 1747 1758 mol->CenterEdge((ofstream *)&cout, &x); 1748 1759 mol->SetBoxDimension(&x); 1760 argptr+=0; 1749 1761 break; 1750 1762 case 'r': … … 1902 1914 if (SaveFlag) 1903 1915 SaveConfig(ConfigFileName, &configuration, periode, molecules); 1904 if ((ExitFlag >= 1)) {1905 delete(mol);1906 delete(periode);1907 return (ExitFlag);1908 }1909 1916 } else { // no arguments, hence scan the elements db 1910 if (periode->LoadPeriodentafel( PathToDatabases))1917 if (periode->LoadPeriodentafel(configuration.databasepath)) 1911 1918 cout << Verbose(0) << "Element list loaded successfully." << endl; 1912 1919 else … … 1914 1921 configuration.RetrieveConfigPathAndName("main_pcp_linux"); 1915 1922 } 1916 return( 0);1923 return(ExitFlag); 1917 1924 }; 1918 1925 … … 1930 1937 ofstream output; 1931 1938 string line; 1932 char ConfigFileName[MAXSTRINGSIZE]; 1933 char *ElementsFileName = NULL; 1939 char *ConfigFileName = NULL; 1934 1940 int j; 1935 1941 1936 1942 // =========================== PARSE COMMAND LINE OPTIONS ==================================== 1937 ConfigFileName[0] = '\0'; 1938 j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName, ElementsFileName); 1939 if (j == 1) return 0; // just for -v and -h options 1940 if (j) return j; // something went wrong 1943 j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName); 1944 switch(j) { 1945 case 0: // something went wrong 1946 delete(molecules); // also free's all molecules contained 1947 delete(periode); 1948 return j; 1949 break; 1950 case 1: // just for -v and -h options 1951 delete(molecules); // also free's all molecules contained 1952 delete(periode); 1953 return 0; 1954 break; 1955 default: 1956 break; 1957 } 1941 1958 1942 1959 // General stuff … … 1952 1969 molecules->insert(mol); 1953 1970 } 1954 if (strlen(ConfigFileName) == 0)1955 strcpy(ConfigFileName, DEFAULTCONFIG);1956 1957 1971 1958 1972 // =========================== START INTERACTIVE SESSION ==================================== … … 2028 2042 2029 2043 // save element data base 2030 if (periode->StorePeriodentafel( ElementsFileName)) //ElementsFileName2044 if (periode->StorePeriodentafel(configuration.databasepath)) //ElementsFileName 2031 2045 cout << Verbose(0) << "Saving of elements.db successful." << endl; 2032 2046 else 2033 2047 cout << Verbose(0) << "Saving of elements.db failed." << endl; 2034 2048 2035 delete(molecules); 2049 delete(molecules); // also free's all molecules contained 2036 2050 delete(periode); 2037 2051 return (0); -
src/config.cpp
r51c910 r36ec71 7 7 #include "molecules.hpp" 8 8 9 /******************************** Functions for class ConfigFileBuffer **********************/ 10 11 /** Structure containing compare function for Ion_Type sorting. 12 */ 13 struct IonTypeCompare { 14 bool operator()(const char* s1, const char *s2) const { 15 char number1[8]; 16 char number2[8]; 17 char *dummy1, *dummy2; 18 //cout << s1 << " " << s2 << endl; 19 dummy1 = strchr(s1, '_')+sizeof(char)*5; // go just after "Ion_Type" 20 dummy2 = strchr(dummy1, '_'); 21 strncpy(number1, dummy1, dummy2-dummy1); // copy the number 22 number1[dummy2-dummy1]='\0'; 23 dummy1 = strchr(s2, '_')+sizeof(char)*5; // go just after "Ion_Type" 24 dummy2 = strchr(dummy1, '_'); 25 strncpy(number2, dummy1, dummy2-dummy1); // copy the number 26 number2[dummy2-dummy1]='\0'; 27 if (atoi(number1) != atoi(number2)) 28 return (atoi(number1) < atoi(number2)); 29 else { 30 dummy1 = strchr(s1, '_')+sizeof(char); 31 dummy1 = strchr(dummy1, '_')+sizeof(char); 32 dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t'); 33 strncpy(number1, dummy1, dummy2-dummy1); // copy the number 34 number1[dummy2-dummy1]='\0'; 35 dummy1 = strchr(s2, '_')+sizeof(char); 36 dummy1 = strchr(dummy1, '_')+sizeof(char); 37 dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t'); 38 strncpy(number2, dummy1, dummy2-dummy1); // copy the number 39 number2[dummy2-dummy1]='\0'; 40 return (atoi(number1) < atoi(number2)); 41 } 42 } 43 }; 44 45 /** Constructor for ConfigFileBuffer class. 46 */ 47 ConfigFileBuffer::ConfigFileBuffer() 48 { 49 NoLines = 0; 50 CurrentLine = 0; 51 buffer = NULL; 52 LineMapping = NULL; 53 } 54 55 /** Constructor for ConfigFileBuffer class with filename to be parsed. 56 * \param *filename file name 57 */ 58 ConfigFileBuffer::ConfigFileBuffer(char *filename) 59 { 60 NoLines = 0; 61 CurrentLine = 0; 62 buffer = NULL; 63 LineMapping = NULL; 64 ifstream *file = NULL; 65 char line[MAXSTRINGSIZE]; 66 67 // prescan number of lines 68 file= new ifstream(filename); 69 if (file == NULL) { 70 cerr << "ERROR: config file " << filename << " missing!" << endl; 71 return; 72 } 73 NoLines = 0; // we're overcounting by one 74 long file_position = file->tellg(); // mark current position 75 do { 76 file->getline(line, 256); 77 NoLines++; 78 } while (!file->eof()); 79 file->clear(); 80 file->seekg(file_position, ios::beg); 81 cout << Verbose(1) << NoLines-1 << " lines were recognized." << endl; 82 83 // allocate buffer's 1st dimension 84 if (buffer != NULL) { 85 cerr << "ERROR: FileBuffer->buffer is not NULL!" << endl; 86 return; 87 } else 88 buffer = (char **) Malloc(sizeof(char *)*NoLines, "ConfigFileBuffer::ConfigFileBuffer: **buffer"); 89 90 // scan each line and put into buffer 91 int lines=0; 92 int i; 93 do { 94 buffer[lines] = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "ConfigFileBuffer::ConfigFileBuffer: *buffer[]"); 95 file->getline(buffer[lines], MAXSTRINGSIZE-1); 96 i = strlen(buffer[lines]); 97 buffer[lines][i] = '\n'; 98 buffer[lines][i+1] = '\0'; 99 lines++; 100 } while((!file->eof()) && (lines < NoLines)); 101 cout << Verbose(1) << lines-1 << " lines were read into the buffer." << endl; 102 103 // close and exit 104 file->close(); 105 file->clear(); 106 delete(file); 107 } 108 109 /** Destructor for ConfigFileBuffer class. 110 */ 111 ConfigFileBuffer::~ConfigFileBuffer() 112 { 113 for(int i=0;i<NoLines;++i) 114 Free((void **)&buffer[i], "ConfigFileBuffer::~ConfigFileBuffer: *buffer[]"); 115 Free((void **)&buffer, "ConfigFileBuffer::~ConfigFileBuffer: **buffer"); 116 Free((void **)&LineMapping, "ConfigFileBuffer::~ConfigFileBuffer: *LineMapping"); 117 } 118 119 120 /** Create trivial mapping. 121 */ 122 void ConfigFileBuffer::InitMapping() 123 { 124 LineMapping = (int *) Malloc(sizeof(int)*NoLines, "ConfigFileBuffer::InitMapping: *LineMapping"); 125 for (int i=0;i<NoLines;i++) 126 LineMapping[i] = i; 127 } 128 129 /** Creates a mapping for the \a *FileBuffer's lines containing the Ion_Type keyword such that they are sorted. 130 * \a *map on return contains a list of NoAtom entries such that going through the list, yields indices to the 131 * lines in \a *FileBuffer in a sorted manner of the Ion_Type?_? keywords. We assume that ConfigFileBuffer::CurrentLine 132 * points to first Ion_Type entry. 133 * \param *FileBuffer pointer to buffer structure 134 * \param NoAtoms of subsequent lines to look at 135 */ 136 void ConfigFileBuffer::MapIonTypesInBuffer(int NoAtoms) 137 { 138 map<const char *, int, IonTypeCompare> LineList; 139 if (LineMapping == NULL) { 140 cerr << "map pointer is NULL: " << LineMapping << endl; 141 return; 142 } 143 144 // put all into hashed map 145 for (int i=0; i<NoAtoms; ++i) { 146 LineList.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i)); 147 } 148 149 // fill map 150 int nr=0; 151 for (map<const char *, int, IonTypeCompare>::iterator runner = LineList.begin(); runner != LineList.end(); ++runner) { 152 if (CurrentLine+nr < NoLines) 153 LineMapping[CurrentLine+(nr++)] = runner->second; 154 else 155 cerr << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl; 156 } 157 } 158 9 159 /************************************* Functions for class config ***************************/ 10 160 … … 14 164 { 15 165 mainname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname"); 16 defaultpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname"); 17 pseudopotpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname"); 18 configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname"); 19 configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname"); 20 ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented"); 21 ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames"); 166 defaultpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: defaultpath"); 167 pseudopotpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: pseudopotpath"); 168 databasepath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: databasepath"); 169 configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: configpath"); 170 configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: configname"); 171 ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "config constructor: *ThermostatImplemented"); 172 ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "config constructor: *ThermostatNames"); 22 173 for (int j=0;j<MaxThermostats;j++) 23 ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), " IonsInitRead: ThermostatNames[]");174 ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "config constructor: ThermostatNames[]"); 24 175 Thermostat = 4; 25 176 alpha = 0.; … … 31 182 configpath[0]='\0'; 32 183 configname[0]='\0'; 33 basis="3-21G"; 34 184 basis = "3-21G"; 35 185 36 186 strcpy(ThermostatNames[0],"None"); … … 116 266 Free((void **)&defaultpath, "config::~config: *defaultpath"); 117 267 Free((void **)&pseudopotpath, "config::~config: *pseudopotpath"); 268 Free((void **)&databasepath, "config::~config: *databasepath"); 118 269 Free((void **)&configpath, "config::~config: *configpath"); 119 270 Free((void **)&configname, "config::~config: *configname"); 271 Free((void **)&ThermostatImplemented, "config::~config: *ThermostatImplemented"); 272 for (int j=0;j<MaxThermostats;j++) 273 Free((void **)&ThermostatNames[j], "config::~config: *ThermostatNames[]"); 274 Free((void **)&ThermostatNames, "config::~config: **ThermostatNames"); 120 275 }; 121 276 … … 550 705 } 551 706 RetrieveConfigPathAndName(filename); 552 // ParseParameters 707 708 // ParseParameterFile 709 struct ConfigFileBuffer *FileBuffer = new ConfigFileBuffer(filename); 710 FileBuffer->InitMapping(); 553 711 554 712 /* Oeffne Hauptparameterdatei */ … … 568 726 /* Namen einlesen */ 569 727 570 ParseForParameter(verbose, file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);571 ParseForParameter(verbose, file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);572 ParseForParameter(verbose, file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);573 ParseForParameter(verbose, file,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);574 ParseForParameter(verbose, file,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);575 576 if (!ParseForParameter(verbose, file,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))728 ParseForParameter(verbose,FileBuffer, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical); 729 ParseForParameter(verbose,FileBuffer, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical); 730 ParseForParameter(verbose,FileBuffer, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical); 731 ParseForParameter(verbose,FileBuffer,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical); 732 ParseForParameter(verbose,FileBuffer,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical); 733 734 if (!ParseForParameter(verbose,FileBuffer,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional)) 577 735 config::Seed = 1; 578 736 579 if(!ParseForParameter(verbose, file,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {737 if(!ParseForParameter(verbose,FileBuffer,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) { 580 738 config::DoOutOrbitals = 0; 581 739 } else { … … 583 741 if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1; 584 742 } 585 ParseForParameter(verbose, file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);743 ParseForParameter(verbose,FileBuffer,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical); 586 744 if (config::DoOutVis < 0) config::DoOutVis = 0; 587 745 if (config::DoOutVis > 1) config::DoOutVis = 1; 588 if (!ParseForParameter(verbose, file,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))746 if (!ParseForParameter(verbose,FileBuffer,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional)) 589 747 config::VectorPlane = -1; 590 if (!ParseForParameter(verbose, file,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))748 if (!ParseForParameter(verbose,FileBuffer,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional)) 591 749 config::VectorCut = 0.; 592 ParseForParameter(verbose, file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);750 ParseForParameter(verbose,FileBuffer,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical); 593 751 if (config::DoOutMes < 0) config::DoOutMes = 0; 594 752 if (config::DoOutMes > 1) config::DoOutMes = 1; 595 if (!ParseForParameter(verbose, file,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))753 if (!ParseForParameter(verbose,FileBuffer,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional)) 596 754 config::DoOutCurrent = 0; 597 755 if (config::DoOutCurrent < 0) config::DoOutCurrent = 0; 598 756 if (config::DoOutCurrent > 1) config::DoOutCurrent = 1; 599 ParseForParameter(verbose, file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);757 ParseForParameter(verbose,FileBuffer,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical); 600 758 if (config::UseAddGramSch < 0) config::UseAddGramSch = 0; 601 759 if (config::UseAddGramSch > 2) config::UseAddGramSch = 2; 602 if(!ParseForParameter(verbose, file,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {760 if(!ParseForParameter(verbose,FileBuffer,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) { 603 761 config::DoWannier = 0; 604 762 } else { … … 606 764 if (config::DoWannier > 1) config::DoWannier = 1; 607 765 } 608 if(!ParseForParameter(verbose, file,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {766 if(!ParseForParameter(verbose,FileBuffer,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) { 609 767 config::CommonWannier = 0; 610 768 } else { … … 612 770 if (config::CommonWannier > 4) config::CommonWannier = 4; 613 771 } 614 if(!ParseForParameter(verbose, file,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {772 if(!ParseForParameter(verbose,FileBuffer,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) { 615 773 config::SawtoothStart = 0.01; 616 774 } else { … … 619 777 } 620 778 621 if (ParseForParameter(verbose, file,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))779 if (ParseForParameter(verbose,FileBuffer,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional)) 622 780 if (config::DoConstrainedMD < 0) 623 781 config::DoConstrainedMD = 0; 624 ParseForParameter(verbose, file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);625 if (!ParseForParameter(verbose, file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))782 ParseForParameter(verbose,FileBuffer,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical); 783 if (!ParseForParameter(verbose,FileBuffer,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional)) 626 784 config::Deltat = 1; 627 ParseForParameter(verbose, file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);628 ParseForParameter(verbose, file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);629 ParseForParameter(verbose, file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);630 //ParseForParameter(verbose, file,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);631 if (!ParseForParameter(verbose, file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))785 ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional); 786 ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional); 787 ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional); 788 //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional); 789 if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional)) 632 790 config::EpsWannier = 1e-8; 633 791 634 792 // stop conditions 635 793 //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1; 636 ParseForParameter(verbose, file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);794 ParseForParameter(verbose,FileBuffer,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical); 637 795 if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3; 638 796 639 ParseForParameter(verbose, file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);640 ParseForParameter(verbose, file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);641 ParseForParameter(verbose, file,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);642 ParseForParameter(verbose, file,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);643 ParseForParameter(verbose, file,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);797 ParseForParameter(verbose,FileBuffer,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical); 798 ParseForParameter(verbose,FileBuffer,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical); 799 ParseForParameter(verbose,FileBuffer,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical); 800 ParseForParameter(verbose,FileBuffer,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical); 801 ParseForParameter(verbose,FileBuffer,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical); 644 802 if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep; 645 803 if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1; 646 804 if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1; 647 805 648 ParseForParameter(verbose, file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);649 ParseForParameter(verbose, file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);650 ParseForParameter(verbose, file,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);651 ParseForParameter(verbose, file,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);652 ParseForParameter(verbose, file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);806 ParseForParameter(verbose,FileBuffer,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical); 807 ParseForParameter(verbose,FileBuffer,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical); 808 ParseForParameter(verbose,FileBuffer,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical); 809 ParseForParameter(verbose,FileBuffer,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical); 810 ParseForParameter(verbose,FileBuffer,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical); 653 811 if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep; 654 812 if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1; … … 656 814 657 815 // Unit cell and magnetic field 658 ParseForParameter(verbose, file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */816 ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */ 659 817 mol->cell_size[0] = BoxLength[0]; 660 818 mol->cell_size[1] = BoxLength[3]; … … 663 821 mol->cell_size[4] = BoxLength[7]; 664 822 mol->cell_size[5] = BoxLength[8]; 665 if (1) fprintf(stderr,"\n");666 667 ParseForParameter(verbose, file,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);668 ParseForParameter(verbose, file,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);669 if (!ParseForParameter(verbose, file,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))823 //if (1) fprintf(stderr,"\n"); 824 825 ParseForParameter(verbose,FileBuffer,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional); 826 ParseForParameter(verbose,FileBuffer,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional); 827 if (!ParseForParameter(verbose,FileBuffer,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional)) 670 828 config::DoFullCurrent = 0; 671 829 if (config::DoFullCurrent < 0) config::DoFullCurrent = 0; … … 678 836 } 679 837 680 ParseForParameter(verbose, file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);681 ParseForParameter(verbose, file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);682 ParseForParameter(verbose, file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);838 ParseForParameter(verbose,FileBuffer,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical); 839 ParseForParameter(verbose,FileBuffer,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical); 840 ParseForParameter(verbose,FileBuffer,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical); 683 841 if (config::Lev0Factor < 2) { 684 842 config::Lev0Factor = 2; 685 843 } 686 ParseForParameter(verbose, file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);844 ParseForParameter(verbose,FileBuffer,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical); 687 845 if (di >= 0 && di < 2) { 688 846 config::RiemannTensor = di; … … 703 861 config::MaxLevel = 3; 704 862 } 705 ParseForParameter(verbose, file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);863 ParseForParameter(verbose,FileBuffer,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical); 706 864 if (config::RiemannLevel < 2) { 707 865 config::RiemannLevel = 2; … … 710 868 config::RiemannLevel = config::MaxLevel-1; 711 869 } 712 ParseForParameter(verbose, file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);870 ParseForParameter(verbose,FileBuffer,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical); 713 871 if (config::LevRFactor < 2) { 714 872 config::LevRFactor = 2; … … 718 876 break; 719 877 } 720 ParseForParameter(verbose, file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);878 ParseForParameter(verbose,FileBuffer,"PsiType", 0, 1, 1, int_type, &di, 1, critical); 721 879 if (di >= 0 && di < 2) { 722 880 config::PsiType = di; … … 727 885 switch (config::PsiType) { 728 886 case 0: // SpinDouble 729 ParseForParameter(verbose, file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);730 ParseForParameter(verbose, file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);887 ParseForParameter(verbose,FileBuffer,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical); 888 ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional); 731 889 break; 732 890 case 1: // SpinUpDown 733 891 if (config::ProcPEGamma % 2) config::ProcPEGamma*=2; 734 ParseForParameter(verbose, file,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);735 ParseForParameter(verbose, file,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);736 ParseForParameter(verbose, file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);892 ParseForParameter(verbose,FileBuffer,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical); 893 ParseForParameter(verbose,FileBuffer,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical); 894 ParseForParameter(verbose,FileBuffer,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional); 737 895 break; 738 896 } … … 740 898 // IonsInitRead 741 899 742 ParseForParameter(verbose, file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);743 ParseForParameter(verbose, file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);744 ParseForParameter(verbose, file,"MaxTypes", 0, 1, 1, int_type, &(config::MaxTypes), 1, critical);745 if (!ParseForParameter(verbose, file,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))900 ParseForParameter(verbose,FileBuffer,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical); 901 ParseForParameter(verbose,FileBuffer,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical); 902 ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(config::MaxTypes), 1, critical); 903 if (!ParseForParameter(verbose,FileBuffer,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional)) 746 904 config::RelativeCoord = 0; 747 if (!ParseForParameter(verbose, file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))905 if (!ParseForParameter(verbose,FileBuffer,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional)) 748 906 config::StructOpt = 0; 749 907 if (MaxTypes == 0) { … … 752 910 // prescan number of ions per type 753 911 cout << Verbose(0) << "Prescanning ions per type: " << endl; 912 int NoAtoms = 0; 754 913 for (int i=0; i < config::MaxTypes; i++) { 755 914 sprintf(name,"Ion_Type%i",i+1); 756 ParseForParameter(verbose, file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);757 ParseForParameter(verbose, file, name, 0, 2, 1, int_type, &Z, 1, critical);915 ParseForParameter(verbose,FileBuffer, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical); 916 ParseForParameter(verbose,FileBuffer, name, 0, 2, 1, int_type, &Z, 1, critical); 758 917 elementhash[i] = periode->FindElement(Z); 759 918 cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 919 NoAtoms += No[i]; 760 920 } 761 921 int repetition = 0; // which repeated keyword shall be read 762 922 923 // sort the lines via the LineMapping 924 sprintf(name,"Ion_Type%i",config::MaxTypes); 925 if (!ParseForParameter(verbose,FileBuffer, (const char*)name, 1, 1, 1, int_type, &value[0], 1, critical)) { 926 cerr << "There are no atoms in the config file!" << endl; 927 return; 928 } 929 FileBuffer->CurrentLine++; 930 //cout << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine]]; 931 FileBuffer->MapIonTypesInBuffer(NoAtoms); 932 //for (int i=0; i<(NoAtoms < 100 ? NoAtoms : 100 < 100 ? NoAtoms : 100);++i) { 933 // cout << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine+i]]; 934 //} 935 763 936 map<int, atom *> AtomList[config::MaxTypes]; 937 map<int, atom *> LinearList; 764 938 if (!FastParsing) { 765 939 // parse in trajectories … … 775 949 neues = new atom(); 776 950 AtomList[i][j] = neues; 951 LinearList[FileBuffer->CurrentLine] = neues; 777 952 neues->type = elementhash[i]; // find element type 778 mol->AddAtom(neues);779 953 } else 780 954 neues = AtomList[i][j]; 781 955 status = (status && 782 ParseForParameter(verbose, file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&783 ParseForParameter(verbose, file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&784 ParseForParameter(verbose, file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&785 ParseForParameter(verbose, file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));956 ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) && 957 ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) && 958 ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) && 959 ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional)); 786 960 if (!status) break; 787 961 … … 799 973 800 974 // parse velocities if present 801 if(!ParseForParameter(verbose, file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))975 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional)) 802 976 neues->v.x[0] = 0.; 803 if(!ParseForParameter(verbose, file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))977 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional)) 804 978 neues->v.x[1] = 0.; 805 if(!ParseForParameter(verbose, file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))979 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional)) 806 980 neues->v.x[2] = 0.; 807 981 for (int d=0;d<NDIM;d++) … … 809 983 810 984 // parse forces if present 811 if(!ParseForParameter(verbose, file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))985 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &value[0], 1,optional)) 812 986 value[0] = 0.; 813 if(!ParseForParameter(verbose, file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))987 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &value[1], 1,optional)) 814 988 value[1] = 0.; 815 if(!ParseForParameter(verbose, file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))989 if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &value[2], 1,optional)) 816 990 value[2] = 0.; 817 991 for (int d=0;d<NDIM;d++) … … 832 1006 repetition++; 833 1007 } 1008 // put atoms into the molecule in their original order 1009 for(map<int, atom*>::iterator runner = LinearList.begin(); runner != LinearList.end(); ++runner) { 1010 mol->AddAtom(runner->second); 1011 } 834 1012 repetition--; 835 1013 cout << "Found " << repetition << " trajectory steps." << endl; … … 838 1016 // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom) 839 1017 repetition = 0; 840 while ( ParseForParameter(verbose, file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&841 ParseForParameter(verbose, file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&842 ParseForParameter(verbose, file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))1018 while ( ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) && 1019 ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) && 1020 ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional)) 843 1021 repetition++; 844 1022 cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl; … … 850 1028 atom *neues = new atom(); 851 1029 // then parse for each atom the coordinates as often as present 852 ParseForParameter(verbose, file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);853 ParseForParameter(verbose, file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);854 ParseForParameter(verbose, file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);855 ParseForParameter(verbose, file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);856 if(!ParseForParameter(verbose, file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))1030 ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical); 1031 ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical); 1032 ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical); 1033 ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical); 1034 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional)) 857 1035 neues->v.x[0] = 0.; 858 if(!ParseForParameter(verbose, file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))1036 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional)) 859 1037 neues->v.x[1] = 0.; 860 if(!ParseForParameter(verbose, file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))1038 if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional)) 861 1039 neues->v.x[2] = 0.; 862 1040 // here we don't care if forces are present (last in trajectories is always equal to current position) … … 867 1045 } 868 1046 } 869 file->close(); 870 delete(file); 1047 delete(FileBuffer); 871 1048 }; 872 1049 … … 1328 1505 * \warning value is modified (both in contents and position)! 1329 1506 * \param verbose 1 - print found value to stderr, 0 - don't 1330 * \param file file to be parsed1507 * \param *file file to be parsed 1331 1508 * \param name Name of value in file (at least 3 chars!) 1332 1509 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning … … 1557 1734 return (found); // true if found, false if not 1558 1735 } 1736 1737 1738 /** Reads parameter from a parsed file. 1739 * The file is either parsed for a certain keyword or if null is given for 1740 * the value in row yth and column xth. If the keyword was necessity#critical, 1741 * then an error is thrown and the programme aborted. 1742 * \warning value is modified (both in contents and position)! 1743 * \param verbose 1 - print found value to stderr, 0 - don't 1744 * \param *FileBuffer pointer to buffer structure 1745 * \param name Name of value in file (at least 3 chars!) 1746 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning 1747 * (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read - 1748 * best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now 1749 * counted from this unresetted position!) 1750 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!) 1751 * \param yth In grid case specifying column number, otherwise the yth \a name matching line 1752 * \param type Type of the Parameter to be read 1753 * \param value address of the value to be read (must have been allocated) 1754 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply 1755 * \param critical necessity of this keyword being specified (optional, critical) 1756 * \return 1 - found, 0 - not found 1757 * \note Routine is taken from the pcp project and hack-a-slack adapted to C++ 1758 */ 1759 int config::ParseForParameter(int verbose, struct ConfigFileBuffer *FileBuffer, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical) { 1760 int i,j; // loop variables 1761 int length = 0, maxlength = -1; 1762 int OldCurrentLine = FileBuffer->CurrentLine; 1763 char *dummy1, *dummy; // pointers in the line that is read in per step 1764 1765 //fprintf(stderr,"Parsing for %s\n",name); 1766 if (repetition == 0) 1767 //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!"); 1768 return 0; 1769 1770 int line = 0; // marks line where parameter was found 1771 int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found 1772 while((found != repetition)) { 1773 dummy1 = dummy = NULL; 1774 do { 1775 dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine++] ]; 1776 if (FileBuffer->CurrentLine >= FileBuffer->NoLines) { 1777 if ((critical) && (found == 0)) { 1778 //Error(InitReading, name); 1779 fprintf(stderr,"Error:InitReading, critical %s not found\n", name); 1780 exit(255); 1781 } else { 1782 FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position 1783 return 0; 1784 } 1785 } 1786 if (dummy1 == NULL) { 1787 if (verbose) fprintf(stderr,"Error reading line %i\n",line); 1788 } else { 1789 //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1); 1790 } 1791 line++; 1792 } while (dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines 1793 1794 // Seek for possible end of keyword on line if given ... 1795 if (name != NULL) { 1796 dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever's nearer 1797 if (dummy == NULL) { 1798 dummy = strchr(dummy1, ' '); // if not found seek for space 1799 while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary 1800 dummy++; 1801 } 1802 if (dummy == NULL) { 1803 dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword) 1804 //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name); 1805 //Free((void **)&free_dummy); 1806 //Error(FileOpenParams, NULL); 1807 } else { 1808 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1); 1809 } 1810 } else dummy = dummy1; 1811 // ... and check if it is the keyword! 1812 //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name)); 1813 if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) { 1814 found++; // found the parameter! 1815 //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy); 1816 1817 if (found == repetition) { 1818 for (i=0;i<xth;i++) { // i = rows 1819 if (type >= grid) { 1820 // grid structure means that grid starts on the next line, not right after keyword 1821 dummy1 = dummy = NULL; 1822 do { 1823 dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[ FileBuffer->CurrentLine++] ]; 1824 if (FileBuffer->CurrentLine >= FileBuffer->NoLines) { 1825 if ((critical) && (found == 0)) { 1826 //Error(InitReading, name); 1827 fprintf(stderr,"Error:InitReading, critical %s not found\n", name); 1828 exit(255); 1829 } else { 1830 FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position 1831 return 0; 1832 } 1833 } 1834 if (dummy1 == NULL) { 1835 if (verbose) fprintf(stderr,"Error reading line %i\n", line); 1836 } else { 1837 //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1); 1838 } 1839 line++; 1840 } while (dummy1 != NULL && (dummy1[0] == '#') || (dummy1[0] == '\n')); 1841 dummy = dummy1; 1842 } else { // simple int, strings or doubles start in the same line 1843 while ((*dummy == '\t') || (*dummy == ' ')) // skip interjacent tabs and spaces 1844 dummy++; 1845 } 1846 1847 for (j=((type >= grid) ? 0 : yth-1);j<yth;j++) { // j = columns 1848 // check for lower triangular area and upper triangular area 1849 if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) { 1850 *((double *)value) = 0.0; 1851 fprintf(stderr,"%f\t",*((double *)value)); 1852 value = (void *)((long)value + sizeof(double)); 1853 //value += sizeof(double); 1854 } else { 1855 // otherwise we must skip all interjacent tabs and spaces and find next value 1856 dummy1 = dummy; 1857 dummy = strchr(dummy1, '\t'); // seek for tab or space 1858 if (dummy == NULL) 1859 dummy = strchr(dummy1, ' '); // if not found seek for space 1860 if (dummy == NULL) { // if still zero returned ... 1861 dummy = strchr(dummy1, '\n'); // ... at line end then 1862 if ((j < yth-1) && (type < 4)) { // check if xth value or not yet 1863 if (critical) { 1864 if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name); 1865 //return 0; 1866 exit(255); 1867 //Error(FileOpenParams, NULL); 1868 } else { 1869 if (!sequential) { // here we need it! 1870 FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position 1871 } 1872 return 0; 1873 } 1874 } 1875 } else { 1876 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy); 1877 } 1878 if (*dummy1 == '#') { 1879 // found comment, skipping rest of line 1880 //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name); 1881 if (!sequential) { // here we need it! 1882 FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position 1883 } 1884 return 0; 1885 } 1886 //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy); 1887 switch(type) { 1888 case (row_int): 1889 *((int *)value) = atoi(dummy1); 1890 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name); 1891 if (verbose) fprintf(stderr,"%i\t",*((int *)value)); 1892 value = (void *)((long)value + sizeof(int)); 1893 //value += sizeof(int); 1894 break; 1895 case(row_double): 1896 case(grid): 1897 case(lower_trigrid): 1898 case(upper_trigrid): 1899 *((double *)value) = atof(dummy1); 1900 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name); 1901 if (verbose) fprintf(stderr,"%lg\t",*((double *)value)); 1902 value = (void *)((long)value + sizeof(double)); 1903 //value += sizeof(double); 1904 break; 1905 case(double_type): 1906 *((double *)value) = atof(dummy1); 1907 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value)); 1908 //value += sizeof(double); 1909 break; 1910 case(int_type): 1911 *((int *)value) = atoi(dummy1); 1912 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value)); 1913 //value += sizeof(int); 1914 break; 1915 default: 1916 case(string_type): 1917 if (value != NULL) { 1918 //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array 1919 maxlength = MAXSTRINGSIZE; 1920 length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum 1921 strncpy((char *)value, dummy1, length); // copy as much 1922 ((char *)value)[length] = '\0'; // and set end marker 1923 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length); 1924 //value += sizeof(char); 1925 } else { 1926 } 1927 break; 1928 } 1929 } 1930 while (*dummy == '\t') 1931 dummy++; 1932 } 1933 } 1934 } 1935 } 1936 } 1937 if ((type >= row_int) && (verbose)) fprintf(stderr,"\n"); 1938 if (!sequential) { 1939 FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position 1940 } 1941 //fprintf(stderr, "End of Parsing\n\n"); 1942 1943 return (found); // true if found, false if not 1944 } -
src/datacreator.cpp
r51c910 r36ec71 54 54 * \return true if file was written successfully 55 55 */ 56 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)56 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 57 57 { 58 58 stringstream filename; … … 87 87 * \return true if file was written successfully 88 88 */ 89 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)89 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 90 90 { 91 91 stringstream filename; … … 124 124 * \return true if file was written successfully 125 125 */ 126 bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateForce)(class MatrixContainer &, int))126 bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix,const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int)) 127 127 { 128 128 stringstream filename; … … 156 156 * \return true if file was written successfully 157 157 */ 158 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateForce)(class MatrixContainer &, int))158 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int)) 159 159 { 160 160 stringstream filename; … … 188 188 * \return true if file was written successfully 189 189 */ 190 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)190 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 191 191 { 192 192 stringstream filename; … … 235 235 * \return true if file was written successfully 236 236 */ 237 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)237 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 238 238 { 239 239 stringstream filename; … … 272 272 * \return true if file was written successfully 273 273 */ 274 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)274 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 275 275 { 276 276 stringstream filename; … … 310 310 * \return true if file was written successfully 311 311 */ 312 bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)312 bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 313 313 { 314 314 stringstream filename; … … 354 354 * \return true if file was written successfully 355 355 */ 356 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum)356 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum) 357 357 { 358 358 stringstream filename; … … 384 384 /** Plot matrix vs. fragment. 385 385 */ 386 bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateFragment)(class MatrixContainer &, int))386 bool CreateDataFragment(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragment)(class MatrixContainer &, int)) 387 387 { 388 388 stringstream filename; … … 449 449 /** Plot matrix vs. fragment. 450 450 */ 451 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int))451 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int)) 452 452 { 453 453 stringstream filename; -
src/datacreator.hpp
r51c910 r36ec71 20 20 bool AppendOutputFile(ofstream &output, const char *dir, const char *filename); 21 21 22 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);23 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);24 bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateForce)(class MatrixContainer &, int));25 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);26 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateForce)(class MatrixContainer &, int));27 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);28 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);29 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);30 bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum);31 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateForce)(class MatrixContainer &, int));32 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, c har *dir, char *prefix, char *msg,char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));22 bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 23 bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 24 bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int)); 25 bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 26 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int)); 27 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 28 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 29 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 30 bool CreateDataDeltaFrobeniusOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum); 31 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateForce)(class MatrixContainer &, int)); 32 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, const char *dir, const char *prefix, const char *msg, const char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int)); 33 33 34 34 void CreateEnergy(class MatrixContainer &Energy, int MatrixNumber); -
src/element.cpp
r51c910 r36ec71 12 12 */ 13 13 element::element() { 14 15 16 17 18 14 Z = -1; 15 No = -1; 16 previous = NULL; 17 next = NULL; 18 sort = NULL; 19 19 }; 20 20 … … 28 28 bool element::Output(ofstream *out) const 29 29 { 30 31 32 //*out << Z << "\t"<< fixed << setprecision(11) << showpoint << mass << "g/mol\t" << name << "\t" << symbol << "\t" << endl;33 34 35 30 if (out != NULL) { 31 *out << name << "\t" << symbol << "\t" << period << "\t" << group << "\t" << block << "\t" << Z << "\t" << mass << "\t" << CovalentRadius << "\t" << VanDerWaalsRadius << endl; 32 //*out << Z << "\t" << fixed << setprecision(11) << showpoint << mass << "g/mol\t" << name << "\t" << symbol << "\t" << endl; 33 return true; 34 } else 35 return false; 36 36 }; 37 37 38 38 /** Prints element data to \a *out. 39 39 * \param *out outstream 40 * \param No 40 * \param No cardinal number of element 41 41 * \param NoOfAtoms total number of atom of this element type 42 42 */ 43 43 bool element::Checkout(ofstream *out, const int Number, const int NoOfAtoms) const 44 44 { 45 46 47 48 49 45 if (out != NULL) { 46 *out << "Ion_Type" << Number << "\t" << NoOfAtoms << "\t" << Z << "\t1.0\t3\t3\t" << fixed << setprecision(11) << showpoint << mass << "\t" << name << "\t" << symbol <<endl; 47 return true; 48 } else 49 return false; 50 50 }; -
src/ellipsoid.cpp
r51c910 r36ec71 2 2 * ellipsoid.cpp 3 3 * 4 * 5 * 4 * Created on: Jan 20, 2009 5 * Author: heber 6 6 */ 7 7 … … 17 17 double SquaredDistanceToEllipsoid(Vector &x, Vector &EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle) 18 18 { 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 19 Vector helper, RefPoint; 20 double distance = -1.; 21 double Matrix[NDIM*NDIM]; 22 double InverseLength[3]; 23 double psi,theta,phi; // euler angles in ZX'Z'' convention 24 25 //cout << Verbose(3) << "Begin of SquaredDistanceToEllipsoid" << endl; 26 27 for(int i=0;i<3;i++) 28 InverseLength[i] = 1./EllipsoidLength[i]; 29 30 // 1. translate coordinate system so that ellipsoid center is in origin 31 helper.CopyVector(&x); 32 helper.SubtractVector(&EllipsoidCenter); 33 RefPoint.CopyVector(&helper); 34 //cout << Verbose(4) << "Translated given point is at " << RefPoint << "." << endl; 35 36 // 2. transform coordinate system by inverse of rotation matrix and of diagonal matrix 37 psi = EllipsoidAngle[0]; 38 theta = EllipsoidAngle[1]; 39 phi = EllipsoidAngle[2]; 40 Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi); 41 Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi); 42 Matrix[2] = sin(psi)*sin(theta); 43 Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi); 44 Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi); 45 Matrix[5] = -cos(psi)*sin(theta); 46 Matrix[6] = sin(theta)*sin(phi); 47 Matrix[7] = sin(theta)*cos(phi); 48 Matrix[8] = cos(theta); 49 helper.MatrixMultiplication(Matrix); 50 helper.Scale(InverseLength); 51 //cout << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl; 52 53 // 3. construct intersection point with unit sphere and ray between origin and x 54 helper.Normalize(); // is simply normalizes vector in distance direction 55 //cout << Verbose(4) << "Transformed intersection is at " << helper << "." << endl; 56 57 // 4. transform back the constructed intersection point 58 psi = -EllipsoidAngle[0]; 59 theta = -EllipsoidAngle[1]; 60 phi = -EllipsoidAngle[2]; 61 helper.Scale(EllipsoidLength); 62 Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi); 63 Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi); 64 Matrix[2] = sin(psi)*sin(theta); 65 Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi); 66 Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi); 67 Matrix[5] = -cos(psi)*sin(theta); 68 Matrix[6] = sin(theta)*sin(phi); 69 Matrix[7] = sin(theta)*cos(phi); 70 Matrix[8] = cos(theta); 71 helper.MatrixMultiplication(Matrix); 72 //cout << Verbose(4) << "Intersection is at " << helper << "." << endl; 73 74 // 5. determine distance between backtransformed point and x 75 distance = RefPoint.DistanceSquared(&helper); 76 //cout << Verbose(4) << "Squared distance between intersection and RefPoint is " << distance << "." << endl; 77 78 return distance; 79 //cout << Verbose(3) << "End of SquaredDistanceToEllipsoid" << endl; 80 80 }; 81 81 … … 83 83 */ 84 84 struct EllipsoidMinimisation { 85 int N;//!< dimension of vector set86 Vector *x;//!< array of vectors85 int N; //!< dimension of vector set 86 Vector *x; //!< array of vectors 87 87 }; 88 88 … … 94 94 double SumSquaredDistance (const gsl_vector * x, void * params) 95 95 { 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 96 Vector *set= ((struct EllipsoidMinimisation *)params)->x; 97 int N = ((struct EllipsoidMinimisation *)params)->N; 98 double SumDistance = 0.; 99 double distance; 100 Vector Center; 101 double EllipsoidLength[3], EllipsoidAngle[3]; 102 103 // put parameters into suitable ellipsoid form 104 for (int i=0;i<3;i++) { 105 Center.x[i] = gsl_vector_get(x, i+0); 106 EllipsoidLength[i] = gsl_vector_get(x, i+3); 107 EllipsoidAngle[i] = gsl_vector_get(x, i+6); 108 } 109 110 // go through all points and sum distance 111 for (int i=0;i<N;i++) { 112 distance = SquaredDistanceToEllipsoid(set[i], Center, EllipsoidLength, EllipsoidAngle); 113 if (!isnan(distance)) { 114 SumDistance += distance; 115 } else { 116 SumDistance = GSL_NAN; 117 break; 118 } 119 } 120 121 //cout << "Current summed distance is " << SumDistance << "." << endl; 122 return SumDistance; 123 123 }; 124 124 … … 132 132 bool FitPointSetToEllipsoid(ofstream *out, Vector *set, int N, Vector *EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle) 133 133 { 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 134 int status = GSL_SUCCESS; 135 *out << Verbose(2) << "Begin of FitPointSetToEllipsoid " << endl; 136 if (N >= 3) { // check that enough points are given (9 d.o.f.) 137 struct EllipsoidMinimisation par; 138 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; 139 gsl_multimin_fminimizer *s = NULL; 140 gsl_vector *ss, *x; 141 gsl_multimin_function minex_func; 142 143 size_t iter = 0; 144 double size; 145 146 /* Starting point */ 147 x = gsl_vector_alloc (9); 148 for (int i=0;i<3;i++) { 149 gsl_vector_set (x, i+0, EllipsoidCenter->x[i]); 150 gsl_vector_set (x, i+3, EllipsoidLength[i]); 151 gsl_vector_set (x, i+6, EllipsoidAngle[i]); 152 } 153 par.x = set; 154 par.N = N; 155 156 /* Set initial step sizes */ 157 ss = gsl_vector_alloc (9); 158 for (int i=0;i<3;i++) { 159 gsl_vector_set (ss, i+0, 0.1); 160 gsl_vector_set (ss, i+3, 1.0); 161 gsl_vector_set (ss, i+6, M_PI/20.); 162 } 163 164 /* Initialize method and iterate */ 165 minex_func.n = 9; 166 minex_func.f = &SumSquaredDistance; 167 minex_func.params = (void *)∥ 168 169 s = gsl_multimin_fminimizer_alloc (T, 9); 170 gsl_multimin_fminimizer_set (s, &minex_func, x, ss); 171 172 do { 173 iter++; 174 status = gsl_multimin_fminimizer_iterate(s); 175 176 if (status) 177 break; 178 179 size = gsl_multimin_fminimizer_size (s); 180 status = gsl_multimin_test_size (size, 1e-2); 181 182 if (status == GSL_SUCCESS) { 183 for (int i=0;i<3;i++) { 184 EllipsoidCenter->x[i] = gsl_vector_get (s->x,i+0); 185 EllipsoidLength[i] = gsl_vector_get (s->x, i+3); 186 EllipsoidAngle[i] = gsl_vector_get (s->x, i+6); 187 } 188 *out << setprecision(3) << Verbose(4) << "Converged fit at: " << *EllipsoidCenter << ", lengths " << EllipsoidLength[0] << ", " << EllipsoidLength[1] << ", " << EllipsoidLength[2] << ", angles " << EllipsoidAngle[0] << ", " << EllipsoidAngle[1] << ", " << EllipsoidAngle[2] << " with summed distance " << s->fval << "." << endl; 189 } 190 191 } while (status == GSL_CONTINUE && iter < 1000); 192 193 gsl_vector_free(x); 194 gsl_vector_free(ss); 195 gsl_multimin_fminimizer_free (s); 196 197 } else { 198 *out << Verbose(3) << "Not enough points provided for fit to ellipsoid." << endl; 199 return false; 200 } 201 *out << Verbose(2) << "End of FitPointSetToEllipsoid" << endl; 202 if (status == GSL_SUCCESS) 203 return true; 204 else 205 return false; 206 206 }; 207 207 … … 213 213 * \param PointsToPick number of points in set to pick 214 214 */ 215 void PickRandomNeighbouredPointSet(ofstream *out, class Tesselation *T, class LinkedCell *LC, Vector *&x, int PointsToPick)216 { 217 int PointsLeft = 0;218 int PointsPicked = 0;219 220 set<int> PickedAtomNrs;// ordered list of picked atoms221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 if (List == NULL) {// set index to it241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 if (PointsLeft < PointsToPick) {// ensure that we can pick enough points in its neighbourhood at all.263 264 265 266 267 268 269 270 current = PickedAtomNrs.find(index);// not present?271 272 273 274 275 } while (PickedAtomNrs.size() < (size_t)PointsToPick);276 277 278 279 280 281 282 283 284 // 285 286 // 287 // 288 // 289 // 290 291 292 293 294 x[PointsPicked++].CopyVector(&(Candidate->x));// we have one more atom picked295 current++;// next pre-picked atom296 297 index++;// next atom nr.298 299 // 300 // 301 302 303 304 for (int i=0;i<PointsPicked;i++)305 306 if (PointsPicked == PointsToPick)// break out of loop if we have all307 308 309 310 215 void PickRandomNeighbouredPointSet(ofstream *out, class Tesselation *T, class LinkedCell *LC, Vector *&x, size_t PointsToPick) 216 { 217 size_t PointsLeft = 0; 218 size_t PointsPicked = 0; 219 int Nlower[NDIM], Nupper[NDIM]; 220 set<int> PickedAtomNrs; // ordered list of picked atoms 221 set<int>::iterator current; 222 int index; 223 atom *Candidate = NULL; 224 LinkedAtoms *List = NULL; 225 *out << Verbose(2) << "Begin of PickRandomPointSet" << endl; 226 227 // allocate array 228 if (x == NULL) { 229 x = new Vector[PointsToPick]; 230 } else { 231 *out << "WARNING: Given pointer to vector array seems already allocated." << endl; 232 } 233 234 do { 235 for(int i=0;i<NDIM;i++) // pick three random indices 236 LC->n[i] = (rand() % LC->N[i]); 237 *out << Verbose(2) << "INFO: Center cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " ... "; 238 // get random cell 239 List = LC->GetCurrentCell(); 240 if (List == NULL) { // set index to it 241 continue; 242 } 243 *out << "with No. " << LC->index << "." << endl; 244 245 *out << Verbose(2) << "LC Intervals:"; 246 for (int i=0;i<NDIM;i++) { 247 Nlower[i] = ((LC->n[i]-1) >= 0) ? LC->n[i]-1 : 0; 248 Nupper[i] = ((LC->n[i]+1) < LC->N[i]) ? LC->n[i]+1 : LC->N[i]-1; 249 *out << " [" << Nlower[i] << "," << Nupper[i] << "] "; 250 } 251 *out << endl; 252 253 // count whether there are sufficient atoms in this cell+neighbors 254 PointsLeft=0; 255 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 256 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 257 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 258 List = LC->GetCurrentCell(); 259 PointsLeft += List->size(); 260 } 261 *out << Verbose(2) << "There are " << PointsLeft << " atoms in this neighbourhood." << endl; 262 if (PointsLeft < PointsToPick) { // ensure that we can pick enough points in its neighbourhood at all. 263 continue; 264 } 265 266 // pre-pick a fixed number of atoms 267 PickedAtomNrs.clear(); 268 do { 269 index = (rand() % PointsLeft); 270 current = PickedAtomNrs.find(index); // not present? 271 if (current == PickedAtomNrs.end()) { 272 //*out << Verbose(2) << "Picking atom nr. " << index << "." << endl; 273 PickedAtomNrs.insert(index); 274 } 275 } while (PickedAtomNrs.size() < PointsToPick); 276 277 index = 0; // now go through all and pick those whose from PickedAtomsNr 278 PointsPicked=0; 279 current = PickedAtomNrs.begin(); 280 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 281 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 282 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 283 List = LC->GetCurrentCell(); 284 // *out << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl; 285 if (List != NULL) { 286 // if (List->begin() != List->end()) 287 // *out << Verbose(2) << "Going through candidates ... " << endl; 288 // else 289 // *out << Verbose(2) << "Cell is empty ... " << endl; 290 for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 291 if ((current != PickedAtomNrs.end()) && (*current == index)) { 292 Candidate = (*Runner); 293 *out << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl; 294 x[PointsPicked++].CopyVector(&(Candidate->x)); // we have one more atom picked 295 current++; // next pre-picked atom 296 } 297 index++; // next atom nr. 298 } 299 // } else { 300 // *out << Verbose(2) << "List for this index not allocated!" << endl; 301 } 302 } 303 *out << Verbose(2) << "The following points were picked: " << endl; 304 for (size_t i=0;i<PointsPicked;i++) 305 *out << Verbose(2) << x[i] << endl; 306 if (PointsPicked == PointsToPick) // break out of loop if we have all 307 break; 308 } while(1); 309 310 *out << Verbose(2) << "End of PickRandomPointSet" << endl; 311 311 }; 312 312 … … 317 317 * \param PointsToPick number of points in set to pick 318 318 */ 319 void PickRandomPointSet(ofstream *out, class Tesselation *T, Vector *&x, int PointsToPick)320 { 321 int PointsLeft =T->PointsOnBoundaryCount;322 int PointsPicked = 0;323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 for (int i=0;i<PointsPicked;i++)350 351 352 319 void PickRandomPointSet(ofstream *out, class Tesselation *T, Vector *&x, size_t PointsToPick) 320 { 321 size_t PointsLeft = (size_t) T->PointsOnBoundaryCount; 322 size_t PointsPicked = 0; 323 double value, threshold; 324 PointMap *List = &T->PointsOnBoundary; 325 *out << Verbose(2) << "Begin of PickRandomPointSet" << endl; 326 327 // allocate array 328 if (x == NULL) { 329 x = new Vector[PointsToPick]; 330 } else { 331 *out << "WARNING: Given pointer to vector array seems already allocated." << endl; 332 } 333 334 if (List != NULL) 335 for (PointMap::iterator Runner = List->begin(); Runner != List->end(); Runner++) { 336 threshold = 1. - (double)(PointsToPick - PointsPicked)/(double)PointsLeft; 337 value = (double)rand()/(double)RAND_MAX; 338 //*out << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": "; 339 if (value > threshold) { 340 x[PointsPicked].CopyVector(&(Runner->second->node->x)); 341 PointsPicked++; 342 //*out << "IN." << endl; 343 } else { 344 //*out << "OUT." << endl; 345 } 346 PointsLeft--; 347 } 348 *out << Verbose(2) << "The following points were picked: " << endl; 349 for (size_t i=0;i<PointsPicked;i++) 350 *out << Verbose(3) << x[i] << endl; 351 352 *out << Verbose(2) << "End of PickRandomPointSet" << endl; 353 353 }; 354 354 … … 363 363 void FindDistributionOfEllipsoids(ofstream *out, class Tesselation *T, class LinkedCell *LCList, int N, int number, const char *filename) 364 364 { 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 EllipsoidCenter.CopyVector(&Center);// use Center of Gravity as initial center of ellipsoid405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 delete[](x);// free allocated memory for point set428 429 430 431 432 433 }; 365 ofstream output; 366 Vector *x = NULL; 367 Vector Center; 368 Vector EllipsoidCenter; 369 double EllipsoidLength[3]; 370 double EllipsoidAngle[3]; 371 double distance, MaxDistance, MinDistance; 372 *out << Verbose(0) << "Begin of FindDistributionOfEllipsoids" << endl; 373 374 // construct center of gravity of boundary point set for initial ellipsoid center 375 Center.Zero(); 376 for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++) 377 Center.AddVector(&Runner->second->node->x); 378 Center.Scale(1./T->PointsOnBoundaryCount); 379 *out << Verbose(1) << "Center is at " << Center << "." << endl; 380 381 // Output header 382 output.open(filename, ios::trunc); 383 output << "# Nr.\tCenterX\tCenterY\tCenterZ\ta\tb\tc\tpsi\ttheta\tphi" << endl; 384 385 // loop over desired number of parameter sets 386 for (;number >0;number--) { 387 *out << Verbose(1) << "Determining data set " << number << " ... " << endl; 388 // pick the point set 389 x = NULL; 390 //PickRandomPointSet(out, T, LCList, x, N); 391 PickRandomNeighbouredPointSet(out, T, LCList, x, N); 392 393 // calculate some sensible starting values for parameter fit 394 MaxDistance = 0.; 395 MinDistance = x[0].ScalarProduct(&x[0]); 396 for (int i=0;i<N;i++) { 397 distance = x[i].ScalarProduct(&x[i]); 398 if (distance > MaxDistance) 399 MaxDistance = distance; 400 if (distance < MinDistance) 401 MinDistance = distance; 402 } 403 //*out << Verbose(2) << "MinDistance " << MinDistance << ", MaxDistance " << MaxDistance << "." << endl; 404 EllipsoidCenter.CopyVector(&Center); // use Center of Gravity as initial center of ellipsoid 405 for (int i=0;i<3;i++) 406 EllipsoidAngle[i] = 0.; 407 EllipsoidLength[0] = sqrt(MaxDistance); 408 EllipsoidLength[1] = sqrt((MaxDistance+MinDistance)/2.); 409 EllipsoidLength[2] = sqrt(MinDistance); 410 411 // fit the parameters 412 if (FitPointSetToEllipsoid(out, x, N, &EllipsoidCenter, &EllipsoidLength[0], &EllipsoidAngle[0])) { 413 *out << Verbose(1) << "Picking succeeded!" << endl; 414 // output obtained parameter set 415 output << number << "\t"; 416 for (int i=0;i<3;i++) 417 output << setprecision(9) << EllipsoidCenter.x[i] << "\t"; 418 for (int i=0;i<3;i++) 419 output << setprecision(9) << EllipsoidLength[i] << "\t"; 420 for (int i=0;i<3;i++) 421 output << setprecision(9) << EllipsoidAngle[i] << "\t"; 422 output << endl; 423 } else { // increase N to pick one more 424 *out << Verbose(1) << "Picking failed!" << endl; 425 number++; 426 } 427 delete[](x); // free allocated memory for point set 428 } 429 // close output and finish 430 output.close(); 431 432 *out << Verbose(0) << "End of FindDistributionOfEllipsoids" << endl; 433 }; -
src/ellipsoid.hpp
r51c910 r36ec71 2 2 * ellipsoid.hpp 3 3 * 4 * 5 * 4 * Created on: Jan 20, 2009 5 * Author: heber 6 6 */ 7 7 … … 22 22 double SumSquaredDistance (const gsl_vector * x, void * params); 23 23 bool FitPointSetToEllipsoid(ofstream *out, Vector *set, int N, Vector *EllipsoidCenter, double *EllipsoidLength, double *EllipsoidAngle); 24 void PickRandomPointSet(ofstream *out, class Tesselation *T, Vector *&x, int N);25 void PickRandomNeighbouredPointSet(ofstream *out, class Tesselation *T, class LinkedCell *LC, Vector *&x, int PointsToPick);24 void PickRandomPointSet(ofstream *out, class Tesselation *T, Vector *&x, size_t PointsToPick); 25 void PickRandomNeighbouredPointSet(ofstream *out, class Tesselation *T, class LinkedCell *LC, Vector *&x, size_t PointsToPick); 26 26 void FindDistributionOfEllipsoids(ofstream *out, class Tesselation *T, class LinkedCell *LCList, int N, int number, const char *filename); 27 27 -
src/helpers.cpp
r51c910 r36ec71 15 15 double ask_value(const char *text) 16 16 { 17 18 19 20 21 22 17 double test = 0.1439851348959832147598734598273456723948652983045928346598365; 18 do { 19 cout << Verbose(0) << text; 20 cin >> test; 21 } while (test == 0.1439851348959832147598734598273456723948652983045928346598365); 22 return test; 23 23 }; 24 24 … … 29 29 #ifdef HAVE_DEBUG 30 30 void debug_in(const char *output, const char *file, const int line) { 31 31 if (output) fprintf(stderr,"DEBUG: in %s at line %i: %s\n", file, line, output); 32 32 } 33 33 #else 34 void debug_in(const char *output, const char *file, const int line) {} 34 void debug_in(const char *output, const char *file, const int line) {} // print nothing 35 35 #endif 36 36 … … 42 42 void * Malloc(size_t size, const char* output) 43 43 { 44 45 46 47 48 44 void *buffer = NULL; 45 buffer = (void *)malloc(size); // alloc 46 if (buffer == NULL) 47 cout << Verbose(0) << "Malloc failed - pointer is NULL: " << output << endl; 48 return(buffer); 49 49 }; 50 50 … … 56 56 void * Calloc(size_t size, const char* output) 57 57 { 58 59 60 61 62 58 void *buffer = NULL; 59 buffer = (void *)calloc(size, (size_t)0); // alloc 60 if (buffer == NULL) 61 cout << Verbose(0) << "Calloc failed - pointer is NULL: " << output << endl; 62 return(buffer); 63 63 }; 64 64 … … 88 88 void Free(void ** buffer, const char* output) 89 89 { 90 91 92 93 94 95 90 if (*buffer == NULL) { 91 //cout << Verbose(5) << "Free not necesary: " << output << endl; 92 } else { 93 free(*buffer); 94 *buffer = NULL; 95 } 96 96 }; 97 97 … … 103 103 char* MallocString(size_t size, const char* output) 104 104 { 105 106 107 108 109 110 for (i=size;i--;)// reset111 112 113 105 size_t i; 106 char *buffer; 107 buffer = (char *)malloc(sizeof(char) * (size+1)); // alloc 108 if (buffer == NULL) 109 cout << Verbose(0) << output << endl; 110 for (i=size;i--;) // reset 111 buffer[i] = i % 2 == 0 ? 'p': 'c'; 112 buffer[size] = '\0'; // and set length marker on its end 113 return(buffer); 114 114 } 115 115 … … 121 121 void bound(double *b, double lower_bound, double upper_bound) 122 122 { 123 124 125 126 127 123 double step = (upper_bound - lower_bound); 124 while (*b >= upper_bound) 125 *b -= step; 126 while (*b < lower_bound) 127 *b += step; 128 128 }; 129 129 … … 134 134 void flip(double *x, double *y) 135 135 { 136 137 138 139 136 double tmp; 137 tmp = *x; 138 *x = *y; 139 *y = tmp; 140 140 }; 141 141 … … 147 147 int pot(int base, int n) 148 148 { 149 150 151 152 153 149 int res = 1; 150 int j; 151 for (j=n;j--;) 152 res *= base; 153 return res; 154 154 }; 155 155 … … 161 161 char *FixedDigitNumber(const int FragmentNumber, const int digits) 162 162 { 163 164 165 166 167 168 169 170 171 172 173 // terminateand fill string array from end backward174 175 176 177 178 179 180 181 163 char *returnstring; 164 int number = FragmentNumber; 165 int order = 0; 166 while (number != 0) { // determine number of digits needed 167 number = (int)floor(((double)number / 10.)); 168 order++; 169 //cout << "Number is " << number << ", order is " << order << "." << endl; 170 } 171 // allocate string 172 returnstring = (char *) Malloc(sizeof(char)*(order+2), "FixedDigitNumber: *returnstring"); 173 // terminate and fill string array from end backward 174 returnstring[order] = '\0'; 175 number = digits; 176 for (int i=order;i--;) { 177 returnstring[i] = '0' + (char)(number % 10); 178 number = (int)floor(((double)number / 10.)); 179 } 180 //cout << returnstring << endl; 181 return returnstring; 182 182 }; 183 183 … … 188 188 bool IsValidNumber( const char *string) 189 189 { 190 191 192 193 194 195 190 int ptr = 0; 191 if ((string[ptr] == '.') || (string[ptr] == '-')) // number may be negative or start with dot 192 ptr++; 193 if ((string[ptr] >= '0') && (string[ptr] <= '9')) 194 return true; 195 return false; 196 196 }; 197 197 -
src/helpers.hpp
r51c910 r36ec71 38 38 39 39 /* Behandelt aufgetretene Fehler. error ist der Fehlertyp(enum Errors) 40 41 40 void *SpecialData ist ein untypisierter Zeiger auf Spezielle Daten zur Fehlerbehandlung. 41 Man koennte auch noch einen Zeiger auf eine Funktion uebergeben */ 42 42 extern void /*@exits@*/ debug(const char *output); 43 43 //__attribute__ ((__return__)); 44 44 #define debug(data) debug_in((data), __FILE__, __LINE__) 45 45 46 46 extern void /*@exits@*/ debug_in(const char *output, 47 48 47 const char *file, const int line); 48 //__attribute__ ((__return__)); 49 49 50 50 double ask_value(const char *text); … … 73 73 template <typename T> bool CreateFatherLookupTable(ofstream *out, T *start, T *end, T **&LookupTable, int count = 0) 74 74 { 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 75 bool status = true; 76 T *Walker = NULL; 77 int AtomNo; 78 79 if (LookupTable != NULL) { 80 *out << "Pointer for Lookup table is not NULL! Aborting ..." <<endl; 81 return false; 82 } 83 84 // count them 85 if (count == 0) { 86 Walker = start; 87 while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron 88 Walker = Walker->next; 89 count = (count < Walker->GetTrueFather()->nr) ? Walker->GetTrueFather()->nr : count; 90 } 91 } 92 if (count <= 0) { 93 *out << "Count of lookup list is 0 or less." << endl; 94 return false; 95 } 96 97 // allocat and fill 98 LookupTable = (T **) Malloc(sizeof(T *)*count, "CreateFatherLookupTable - **LookupTable"); 99 if (LookupTable == NULL) { 100 cerr << "LookupTable memory allocation failed!" << endl; 101 status = false; 102 } else { 103 for (int i=0;i<count;i++) 104 LookupTable[i] = NULL; 105 Walker = start; 106 while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron 107 Walker = Walker->next; 108 AtomNo = Walker->GetTrueFather()->nr; 109 if ((AtomNo >= 0) && (AtomNo < count)) { 110 //*out << "Setting LookupTable[" << AtomNo << "] to " << *Walker << endl; 111 LookupTable[AtomNo] = Walker; 112 } else { 113 *out << "Walker " << *Walker << " exceeded range of nuclear ids [0, " << count << ")." << endl; 114 status = false; 115 break; 116 } 117 } 118 } 119 120 return status; 121 121 }; 122 122 … … 129 129 template <typename X> void link(X *walker, X *end) 130 130 { 131 132 133 134 135 136 131 X *vorher = end->previous; 132 if (vorher != NULL) 133 vorher->next = walker; 134 end->previous = walker; 135 walker->previous = vorher; 136 walker->next = end; 137 137 }; 138 138 … … 143 143 template <typename X> void unlink(X *walker) 144 144 { 145 146 147 148 145 if (walker->next != NULL) 146 walker->next->previous = walker->previous; 147 if (walker->previous != NULL) 148 walker->previous->next = walker->next; 149 149 }; 150 150 151 151 /** Adds new item before an item \a *end in a list. 152 * \param *pointer 153 * \param *end 152 * \param *pointer item to be added 153 * \param *end end of list 154 154 * \return true - addition succeeded, false - unable to add item to list 155 155 */ 156 template <typename X> 157 { 158 159 160 161 162 163 164 156 template <typename X> bool add(X *pointer, X *end) 157 { 158 if (end != NULL) { 159 link(pointer, end); 160 } else { 161 pointer->previous = NULL; 162 pointer->next = NULL; 163 } 164 return true; 165 165 }; 166 166 167 167 /** Finds item in list 168 * \param *suche 169 * \param *start 170 * \param *end 168 * \param *suche search criteria 169 * \param *start begin of list 170 * \param *end end of list 171 171 * \return X - if found, NULL - if not found 172 172 */ 173 173 template <typename X, typename Y> X * find(Y *suche, X *start, X *end) 174 174 { 175 176 177 178 179 180 175 X *walker = start; 176 while (walker->next != end) { // go through list 177 walker = walker->next; // step onward beforehand 178 if (*walker->sort == *suche) return (walker); 179 } 180 return NULL; 181 181 }; 182 182 … … 187 187 template <typename X> void removewithoutcheck(X *walker) 188 188 { 189 190 191 192 193 189 if (walker != NULL) { 190 unlink(walker); 191 delete(walker); 192 walker = NULL; 193 } 194 194 }; 195 195 196 196 /** Removes an item from the list, checks if exists. 197 197 * Checks beforehand if atom is really within molecule list. 198 * \param *pointer 199 * \param *start 200 * \param *end 198 * \param *pointer item to be removed 199 * \param *start begin of list 200 * \param *end end of list 201 201 * \return true - removing succeeded, false - given item not found in list 202 202 */ 203 203 template <typename X> bool remove(X *pointer, X *start, X *end) 204 204 { 205 206 /* 207 208 if (walker == end) return false;// item not found in list209 210 211 212 213 214 215 205 X *walker = find (pointer->sort, start, end); 206 /* while (walker->next != pointer) { // search through list 207 walker = walker->next; 208 if (walker == end) return false; // item not found in list 209 }*/ 210 // atom found, now unlink 211 if (walker != NULL) 212 removewithoutcheck(walker); 213 else 214 return false; 215 return true; 216 216 }; 217 217 … … 223 223 template <typename X> bool cleanup(X *start, X *end) 224 224 { 225 226 227 228 229 230 231 232 233 234 235 225 X *pointer = start->next; 226 X *walker; 227 while (pointer != end) { // go through list 228 walker = pointer; // mark current 229 pointer = pointer->next; // step onward beforehand 230 // remove walker 231 unlink(walker); 232 delete(walker); 233 walker = NULL; 234 } 235 return true; 236 236 }; 237 237 … … 242 242 template <typename X> X *GetFirst(X *me) 243 243 { 244 245 246 247 244 X *Binder = me; 245 while(Binder->previous != NULL) 246 Binder = Binder->previous; 247 return Binder; 248 248 }; 249 249 … … 254 254 template <typename X> X *GetLast(X *me) 255 255 { 256 257 258 259 256 X *Binder = me; 257 while(Binder->next != NULL) 258 Binder = Binder->next; 259 return Binder; 260 260 }; 261 261 … … 266 266 template <typename X> void Free2DArray(X **ptr, int dim) 267 267 { 268 269 270 271 272 273 274 268 int i; 269 if (ptr != NULL) { 270 for(i=dim;i--;) 271 if (ptr[i] != NULL) 272 free(ptr[i]); 273 free(ptr); 274 } 275 275 }; 276 276 … … 281 281 class Verbose 282 282 { 283 284 285 286 287 288 283 public: 284 Verbose(int value) : Verbosity(value) { } 285 286 ostream& print (ostream &ost) const; 287 private: 288 int Verbosity; 289 289 }; 290 290 … … 295 295 class Binary 296 296 { 297 298 299 300 301 302 297 public: 298 Binary(int value) : BinaryNumber(value) { } 299 300 ostream& print (ostream &ost) const; 301 private: 302 int BinaryNumber; 303 303 }; 304 304 -
src/joiner.cpp
r51c910 r36ec71 110 110 // ---------- Parse the KeySets into an array --------------- 111 111 if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1; 112 112 113 if (!KeySet.ParseManyBodyTerms()) return 1; 113 114 -
src/linkedcell.cpp
r51c910 r36ec71 6 6 LinkedCell::LinkedCell() 7 7 { 8 9 10 11 12 13 14 8 LC = NULL; 9 for(int i=0;i<NDIM;i++) 10 N[i] = 0; 11 index = -1; 12 RADIUS = 0.; 13 max.Zero(); 14 min.Zero(); 15 15 }; 16 16 … … 21 21 LinkedCell::LinkedCell(molecule *mol, double radius) 22 22 { 23 23 atom *Walker = NULL; 24 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 25 RADIUS = radius; 26 LC = NULL; 27 for(int i=0;i<NDIM;i++) 28 N[i] = 0; 29 index = -1; 30 max.Zero(); 31 min.Zero(); 32 cout << Verbose(1) << "Begin of LinkedCell" << endl; 33 if (mol->start->next == mol->end) { 34 cerr << "ERROR: molecule contains no atoms!" << endl; 35 return; 36 } 37 // 1. find max and min per axis of atoms 38 Walker = mol->start->next; 39 for (int i=0;i<NDIM;i++) { 40 max.x[i] = Walker->x.x[i]; 41 min.x[i] = Walker->x.x[i]; 42 } 43 while (Walker != mol->end) { 44 for (int i=0;i<NDIM;i++) { 45 if (max.x[i] < Walker->x.x[i]) 46 max.x[i] = Walker->x.x[i]; 47 if (min.x[i] > Walker->x.x[i]) 48 min.x[i] = Walker->x.x[i]; 49 } 50 Walker = Walker->next; 51 } 52 cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl; 53 53 54 55 56 57 58 54 // 2. find then umber of cells per axis 55 for (int i=0;i<NDIM;i++) { 56 N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1; 57 } 58 cout << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl; 59 59 60 // 3. allocate the lists 61 cout << Verbose(2) << "Allocating cells ... " << endl; 62 if (LC != NULL) { 63 cout << Verbose(1) << "ERROR: Linked Cell list is already allocated, I do nothing." << endl; 64 return; 65 } 66 LC = new LinkedAtoms[N[0]*N[1]*N[2]]; 67 for (index=0;index<N[0]*N[1]*N[2];index++) { 68 LC [index].clear(); 69 } 60 // 3. allocate the lists 61 cout << Verbose(2) << "Allocating cells ... "; 62 if (LC != NULL) { 63 cout << Verbose(1) << "ERROR: Linked Cell list is already allocated, I do nothing." << endl; 64 return; 65 } 66 LC = new LinkedAtoms[N[0]*N[1]*N[2]]; 67 for (index=0;index<N[0]*N[1]*N[2];index++) { 68 LC [index].clear(); 69 } 70 cout << "done." << endl; 70 71 71 // 4. put each atom into its respective cell 72 Walker = mol->start; 73 while (Walker->next != mol->end) { 74 Walker = Walker->next; 75 for (int i=0;i<NDIM;i++) { 76 n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS); 77 } 78 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2]; 79 LC[index].push_back(Walker); 80 cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl; 81 } 82 cout << Verbose(1) << "End of LinkedCell" << endl; 72 // 4. put each atom into its respective cell 73 cout << Verbose(2) << "Filling cells ... "; 74 Walker = mol->start; 75 while (Walker->next != mol->end) { 76 Walker = Walker->next; 77 for (int i=0;i<NDIM;i++) { 78 n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS); 79 } 80 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2]; 81 LC[index].push_back(Walker); 82 //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl; 83 } 84 cout << "done." << endl; 85 cout << Verbose(1) << "End of LinkedCell" << endl; 83 86 }; 84 87 … … 87 90 LinkedCell::~LinkedCell() 88 91 { 89 90 91 92 93 94 95 96 97 92 if (LC != NULL) 93 for (index=0;index<N[0]*N[1]*N[2];index++) 94 LC[index].clear(); 95 delete[](LC); 96 for(int i=0;i<NDIM;i++) 97 N[i] = 0; 98 index = -1; 99 max.Zero(); 100 min.Zero(); 98 101 }; 99 102 … … 103 106 bool LinkedCell::CheckBounds() 104 107 { 105 106 107 108 109 110 108 bool status = true; 109 for(int i=0;i<NDIM;i++) 110 status = status && ((n[i] >=0) && (n[i] < N[i])); 111 if (!status) 112 cerr << "ERROR: indices are out of bounds!" << endl; 113 return status; 111 114 }; 112 115 … … 117 120 LinkedAtoms* LinkedCell::GetCurrentCell() 118 121 { 119 120 121 122 123 124 122 if (CheckBounds()) { 123 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2]; 124 return (&(LC[index])); 125 } else { 126 return NULL; 127 } 125 128 }; 126 129 … … 131 134 bool LinkedCell::SetIndexToAtom(atom *Walker) 132 135 { 133 134 135 136 137 138 139 140 141 142 143 144 145 136 bool status = false; 137 for (int i=0;i<NDIM;i++) { 138 n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS); 139 } 140 index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2]; 141 if (CheckBounds()) { 142 for (LinkedAtoms::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++) 143 status = status || ((*Runner) == Walker); 144 return status; 145 } else { 146 cerr << Verbose(1) << "ERROR: Atom "<< *Walker << " at " << Walker->x << " is out of bounds." << endl; 147 return false; 148 } 146 149 }; 147 150 … … 152 155 bool LinkedCell::SetIndexToVector(Vector *x) 153 156 { 154 155 156 157 158 159 160 161 162 157 bool status = true; 158 for (int i=0;i<NDIM;i++) { 159 n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS); 160 if (max.x[i] < x->x[i]) 161 status = false; 162 if (min.x[i] > x->x[i]) 163 status = false; 164 } 165 return status; 163 166 }; 164 167 -
src/linkedcell.hpp
r51c910 r36ec71 12 12 13 13 class LinkedCell{ 14 15 Vector max;// upper boundary16 Vector min;// lower boundary17 LinkedAtoms *LC;// linked cell list18 double RADIUS;// cell edge length19 int N[NDIM];// number of cells per axis20 int n[NDIM];// temporary variable for current cell per axis21 int index;// temporary index variable , access by index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];14 public: 15 Vector max; // upper boundary 16 Vector min; // lower boundary 17 LinkedAtoms *LC; // linked cell list 18 double RADIUS; // cell edge length 19 int N[NDIM]; // number of cells per axis 20 int n[NDIM]; // temporary variable for current cell per axis 21 int index; // temporary index variable , access by index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2]; 22 22 23 24 25 26 27 28 29 23 LinkedCell(); 24 LinkedCell(molecule *mol, double RADIUS); 25 ~LinkedCell(); 26 LinkedAtoms* GetCurrentCell(); 27 bool SetIndexToAtom(atom *Walker); 28 bool SetIndexToVector(Vector *x); 29 bool CheckBounds(); 30 30 31 32 33 34 31 // not implemented yet 32 bool AddAtom(atom *Walker); 33 bool DeleteAtom(atom *Walker); 34 bool MoveAtom(atom *Walker); 35 35 }; 36 36 -
src/molecules.cpp
r51c910 r36ec71 632 632 /** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths. 633 633 * \param *out output stream for debugging 634 * \param *BoxLengths box lengths 635 */ 636 bool molecule::CenterInBox(ofstream *out, Vector *BoxLengths) 634 */ 635 bool molecule::CenterInBox(ofstream *out) 637 636 { 638 637 bool status = true; 639 638 atom *ptr = NULL; 640 Vector *min = new Vector; 641 Vector *max = new Vector; 642 643 // gather min and max for each axis 639 Vector x; 640 double *M = ReturnFullMatrixforSymmetric(cell_size); 641 double *Minv = x.InverseMatrix(M); 642 double value; 643 644 // cout << "The box matrix is :" << endl; 645 // for (int i=0;i<NDIM;++i) { 646 // for (int j=0;j<NDIM;++j) 647 // cout << M[i*NDIM+j] << "\t"; 648 // cout << endl; 649 // } 650 // cout << "And its inverse is :" << endl; 651 // for (int i=0;i<NDIM;++i) { 652 // for (int j=0;j<NDIM;++j) 653 // cout << Minv[i*NDIM+j] << "\t"; 654 // cout << endl; 655 // } 656 // go through all atoms 644 657 ptr = start->next; // start at first in list 645 if (ptr != end) { //list not empty? 646 for (int i=NDIM;i--;) { 647 max->x[i] = ptr->x.x[i]; 648 min->x[i] = ptr->x.x[i]; 649 } 658 if (ptr != end) { //list not empty? 650 659 while (ptr->next != end) { // continue with second if present 651 660 ptr = ptr->next; 652 661 //ptr->Output(1,1,out); 653 for (int i=NDIM;i--;) { 654 max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i]; 655 min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i]; 656 } 657 } 658 } 659 // sanity check 660 for(int i=NDIM;i--;) { 661 if (max->x[i] - min->x[i] > BoxLengths->x[i]) 662 status = false; 663 } 664 // warn if check failed 665 if (!status) 666 *out << "WARNING: molecule is bigger than defined box!" << endl; 667 else { // else center in box 668 max->AddVector(min); 669 max->Scale(-1.); 670 max->AddVector(BoxLengths); 671 max->Scale(0.5); 672 Translate(max); 673 Center.Zero(); 674 } 675 676 // free and exit 677 delete(min); 678 delete(max); 662 // multiply its vector with matrix inverse 663 x.CopyVector(&ptr->x); 664 x.MatrixMultiplication(Minv); 665 // truncate to [0,1] for each axis 666 for (int i=0;i<NDIM;i++) { 667 value = floor(x.x[i]); // next lower integer 668 if (x.x[i] >=0) { 669 x.x[i] -= value; 670 } else { 671 x.x[i] += value+1; 672 } 673 } 674 } 675 } 676 delete(M); 677 delete(Minv); 679 678 return status; 680 679 }; … … 806 805 }; 807 806 807 808 808 /** Centers the center of gravity of the atoms at (0,0,0). 809 809 * \param *out output stream for debugging … … 814 814 Center.CopyVector(newcenter); 815 815 }; 816 816 817 817 818 /** Scales all atoms by \a *factor. … … 844 845 } 845 846 }; 847 848 /** Translate the molecule periodically in the box. 849 * \param trans[] translation vector. 850 */ 851 void molecule::TranslatePeriodically(const Vector *trans) 852 { 853 atom *ptr = NULL; 854 Vector x; 855 double *M = ReturnFullMatrixforSymmetric(cell_size); 856 double *Minv = x.InverseMatrix(M); 857 double value; 858 859 // go through all atoms 860 ptr = start->next; // start at first in list 861 while (ptr != end) { // continue with second if present 862 //ptr->Output(1,1,out); 863 // multiply its vector with matrix inverse 864 x.CopyVector(&ptr->x); 865 x.Translate(trans); 866 x.MatrixMultiplication(Minv); 867 // truncate to [0,1] for each axis 868 for (int i=0;i<NDIM;i++) { 869 value = floor(fabs(x.x[i])); // next lower integer 870 if (x.x[i] >=0) { 871 x.x[i] -= value; 872 } else { 873 x.x[i] += value+1; 874 } 875 } 876 // matrix multiply 877 x.MatrixMultiplication(M); 878 ptr->x.CopyVector(&x); 879 for (int j=0;j<MDSteps;j++) { 880 x.CopyVector(&Trajectories[ptr].R.at(j)); 881 x.Translate(trans); 882 x.MatrixMultiplication(Minv); 883 // truncate to [0,1] for each axis 884 for (int i=0;i<NDIM;i++) { 885 value = floor(x.x[i]); // next lower integer 886 if (x.x[i] >=0) { 887 x.x[i] -= value; 888 } else { 889 x.x[i] += value+1; 890 } 891 } 892 // matrix multiply 893 x.MatrixMultiplication(M); 894 Trajectories[ptr].R.at(j).CopyVector(&x); 895 } 896 ptr = ptr->next; 897 } 898 delete(M); 899 delete(Minv); 900 }; 901 846 902 847 903 /** Mirrors all atoms against a given plane. … … 873 929 874 930 do { 875 center.Zero();931 Center.Zero(); 876 932 flag = true; 877 933 while (Walker->next != end) { … … 900 956 Testvector.AddVector(&Translationvector); 901 957 Testvector.MatrixMultiplication(matrix); 902 center.AddVector(&Testvector);958 Center.AddVector(&Testvector); 903 959 cout << Verbose(1) << "vector is: "; 904 960 Testvector.Output((ofstream *)&cout); … … 913 969 Testvector.AddVector(&Translationvector); 914 970 Testvector.MatrixMultiplication(matrix); 915 center.AddVector(&Testvector);971 Center.AddVector(&Testvector); 916 972 cout << Verbose(1) << "Hydrogen vector is: "; 917 973 Testvector.Output((ofstream *)&cout); … … 924 980 } while (!flag); 925 981 Free((void **)&matrix, "molecule::DetermineCenter: *matrix"); 926 center.Scale(1./(double)AtomCount);982 Center.Scale(1./(double)AtomCount); 927 983 }; 928 984 … … 937 993 Vector *CenterOfGravity = DetermineCenterOfGravity(out); 938 994 939 Center AtVector(out, CenterOfGravity);995 CenterPeriodic(out); 940 996 941 997 // reset inertia tensor … … 2061 2117 bool molecule::Output(ofstream *out) 2062 2118 { 2063 element *runner;2064 2119 atom *walker = NULL; 2065 int ElementNo , AtomNo;2120 int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS]; 2066 2121 CountElements(); 2067 2122 2123 for (int i=0;i<MAX_ELEMENTS;++i) { 2124 AtomNo[i] = 0; 2125 ElementNo[i] = 0; 2126 } 2068 2127 if (out == NULL) { 2069 2128 return false; 2070 2129 } else { 2071 2130 *out << "#Ion_TypeNr._Nr.R[0] R[1] R[2] MoveType (0 MoveIon, 1 FixedIon)" << endl; 2072 ElementNo = 0; 2073 runner = elemente->start; 2074 while (runner->next != elemente->end) { // go through every element 2075 runner = runner->next; 2076 if (ElementsInMolecule[runner->Z]) { // if this element got atoms 2077 ElementNo++; 2078 AtomNo = 0; 2079 walker = start; 2080 while (walker->next != end) { // go through every atom of this element 2081 walker = walker->next; 2082 if (walker->type == runner) { // if this atom fits to element 2083 AtomNo++; 2084 walker->Output(ElementNo, AtomNo, out); // removed due to trajectories 2085 } 2086 } 2087 } 2131 walker = start; 2132 while (walker->next != end) { // go through every atom of this element 2133 walker = walker->next; 2134 ElementNo[walker->type->Z] = 1; 2135 } 2136 int current=1; 2137 for (int i=0;i<MAX_ELEMENTS;++i) { 2138 if (ElementNo[i] == 1) 2139 ElementNo[i] = current++; 2140 } 2141 walker = start; 2142 while (walker->next != end) { // go through every atom of this element 2143 walker = walker->next; 2144 AtomNo[walker->type->Z]++; 2145 walker->Output(ElementNo[walker->type->Z], AtomNo[walker->type->Z], out); // removed due to trajectories 2088 2146 } 2089 2147 return true; … … 2096 2154 bool molecule::OutputTrajectories(ofstream *out) 2097 2155 { 2098 element *runner = NULL;2099 2156 atom *walker = NULL; 2100 int ElementNo , AtomNo;2157 int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS]; 2101 2158 CountElements(); 2102 2159 … … 2110 2167 *out << "# ====== MD step " << step << " =========" << endl; 2111 2168 } 2112 ElementNo = 0; 2113 runner = elemente->start; 2114 while (runner->next != elemente->end) { // go through every element 2115 runner = runner->next; 2116 if (ElementsInMolecule[runner->Z]) { // if this element got atoms 2117 ElementNo++; 2118 AtomNo = 0; 2119 walker = start; 2120 while (walker->next != end) { // go through every atom of this element 2121 walker = walker->next; 2122 if (walker->type == runner) { // if this atom fits to element 2123 AtomNo++; 2124 *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t" << fixed << setprecision(9) << showpoint; 2125 *out << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2]; 2126 *out << "\t" << walker->FixedIon; 2127 if (Trajectories[walker].U.at(step).Norm() > MYEPSILON) 2128 *out << "\t" << scientific << setprecision(6) << Trajectories[walker].U.at(step).x[0] << "\t" << Trajectories[walker].U.at(step).x[1] << "\t" << Trajectories[walker].U.at(step).x[2] << "\t"; 2129 if (Trajectories[walker].F.at(step).Norm() > MYEPSILON) 2130 *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t"; 2131 *out << "\t# Number in molecule " << walker->nr << endl; 2132 } 2133 } 2134 } 2169 for (int i=0;i<MAX_ELEMENTS;++i) { 2170 AtomNo[i] = 0; 2171 ElementNo[i] = 0; 2172 } 2173 walker = start; 2174 while (walker->next != end) { // go through every atom of this element 2175 walker = walker->next; 2176 ElementNo[walker->type->Z] = 1; 2177 } 2178 int current=1; 2179 for (int i=0;i<MAX_ELEMENTS;++i) { 2180 if (ElementNo[i] == 1) 2181 ElementNo[i] = current++; 2182 } 2183 walker = start; 2184 while (walker->next != end) { // go through every atom of this element 2185 walker = walker->next; 2186 AtomNo[walker->type->Z]++; 2187 *out << "Ion_Type" << ElementNo[walker->type->Z] << "_" << AtomNo[walker->type->Z] << "\t" << fixed << setprecision(9) << showpoint; 2188 *out << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2]; 2189 *out << "\t" << walker->FixedIon; 2190 if (Trajectories[walker].U.at(step).Norm() > MYEPSILON) 2191 *out << "\t" << scientific << setprecision(6) << Trajectories[walker].U.at(step).x[0] << "\t" << Trajectories[walker].U.at(step).x[1] << "\t" << Trajectories[walker].U.at(step).x[2] << "\t"; 2192 if (Trajectories[walker].F.at(step).Norm() > MYEPSILON) 2193 *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t"; 2194 *out << "\t# Number in molecule " << walker->nr << endl; 2135 2195 } 2136 2196 } … … 2205 2265 { 2206 2266 atom *walker = NULL; 2207 int AtomNo = 0 , ElementNo;2267 int AtomNo = 0; 2208 2268 time_t now; 2209 element *runner = NULL;2210 2269 2211 2270 now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' … … 2217 2276 if (out != NULL) { 2218 2277 *out << AtomNo << "\n\tCreated by molecuilder on " << ctime(&now); 2219 ElementNo = 0; 2220 runner = elemente->start; 2221 while (runner->next != elemente->end) { // go through every element 2222 runner = runner->next; 2223 if (ElementsInMolecule[runner->Z]) { // if this element got atoms 2224 ElementNo++; 2225 walker = start; 2226 while (walker->next != end) { // go through every atom of this element 2227 walker = walker->next; 2228 if (walker->type == runner) { // if this atom fits to element 2229 walker->OutputXYZLine(out); 2230 } 2231 } 2232 } 2278 walker = start; 2279 while (walker->next != end) { // go through every atom of this element 2280 walker = walker->next; 2281 walker->OutputXYZLine(out); 2233 2282 } 2234 2283 return true; … … 2451 2500 Vector x; 2452 2501 int FalseBondDegree = 0; 2453 2502 2454 2503 BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem); 2455 2504 *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl; … … 4946 4995 matrix[0] = symm[0]; 4947 4996 matrix[1] = symm[1]; 4948 matrix[2] = symm[ 3];4997 matrix[2] = symm[2]; 4949 4998 matrix[3] = symm[1]; 4950 matrix[4] = symm[ 2];4999 matrix[4] = symm[3]; 4951 5000 matrix[5] = symm[4]; 4952 matrix[6] = symm[ 3];5001 matrix[6] = symm[2]; 4953 5002 matrix[7] = symm[4]; 4954 5003 matrix[8] = symm[5]; -
src/molecules.hpp
r51c910 r36ec71 59 59 #define PointPair pair < int, class BoundaryPointSet * > 60 60 #define PointTestPair pair < PointMap::iterator, bool > 61 62 #define LineMap map < int, class BoundaryLineSet * > 61 #define CandidateList list <class CandidateForTesselation *> 62 63 #define LineMap multimap < int, class BoundaryLineSet * > 63 64 #define LinePair pair < int, class BoundaryLineSet * > 64 65 #define LineTestPair pair < LineMap::iterator, bool > … … 73 74 #define MoleculeList list <molecule *> 74 75 #define MoleculeListTest pair <MoleculeList::iterator, bool> 76 77 #define LinkedAtoms list <atom *> 75 78 76 79 /******************************** Some small functions and/or structures **********************************/ … … 257 260 void CountElements(); 258 261 void CalculateOrbitals(class config &configuration); 259 bool CenterInBox(ofstream *out , Vector *BoxLengths);262 bool CenterInBox(ofstream *out); 260 263 void CenterEdge(ofstream *out, Vector *max); 261 264 void CenterOrigin(ofstream *out); … … 263 266 void CenterAtVector(ofstream *out, Vector *newcenter); 264 267 void Translate(const Vector *x); 268 void TranslatePeriodically(const Vector *trans); 265 269 void Mirror(const Vector *x); 266 270 void Align(Vector *n); … … 403 407 }; 404 408 409 class ConfigFileBuffer { 410 public: 411 char **buffer; 412 int *LineMapping; 413 int CurrentLine; 414 int NoLines; 415 416 ConfigFileBuffer(); 417 ConfigFileBuffer(char *filename); 418 ~ConfigFileBuffer(); 419 420 void InitMapping(); 421 void MapIonTypesInBuffer(int NoAtoms); 422 }; 423 424 405 425 /** The config file. 406 426 * The class contains all parameters that control a dft run also functions to load and save. … … 422 442 string basis; 423 443 444 char *databasepath; 445 424 446 int DoConstrainedMD; 425 447 int MaxOuterStep; … … 433 455 int ScaleTempStep; 434 456 435 457 private: 436 458 char *mainname; 437 459 char *defaultpath; … … 486 508 487 509 int ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical); 510 int ParseForParameter(int verbose, struct ConfigFileBuffer *FileBuffer, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical); 488 511 489 512 public: -
src/parser.cpp
r51c910 r36ec71 249 249 * \return parsing successful 250 250 */ 251 bool MatrixContainer::ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns)251 bool MatrixContainer::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns) 252 252 { 253 253 char filename[1023]; … … 562 562 * \return parsing successful 563 563 */ 564 bool EnergyMatrix::ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns)564 bool EnergyMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns) 565 565 { 566 566 char filename[1024]; … … 598 598 * \return parsing successful 599 599 */ 600 bool ForceMatrix::ParseIndices(c har *name)600 bool ForceMatrix::ParseIndices(const char *name) 601 601 { 602 602 ifstream input; … … 676 676 * \return parsing successful 677 677 */ 678 bool ForceMatrix::ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns)678 bool ForceMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns) 679 679 { 680 680 char filename[1023]; … … 902 902 * \return parsing successful 903 903 */ 904 bool HessianMatrix::ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns)904 bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns) 905 905 { 906 906 char filename[1023]; -
src/parser.hpp
r51c910 r36ec71 59 59 bool InitialiseIndices(class MatrixContainer *Matrix = NULL); 60 60 bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr); 61 virtual bool ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns);61 virtual bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns); 62 62 bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter); 63 63 bool ResetMatrix(); … … 79 79 bool ParseIndices(); 80 80 bool SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign); 81 bool ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns);81 bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns); 82 82 }; 83 83 … … 86 86 class ForceMatrix : public MatrixContainer { 87 87 public: 88 bool ParseIndices(c har *name);88 bool ParseIndices(const char *name); 89 89 bool SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign); 90 bool ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns);90 bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns); 91 91 }; 92 92 … … 100 100 bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order); 101 101 bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign); 102 bool ParseFragmentMatrix(c har *name,char *prefix, string suffix, int skiplines, int skipcolumns);102 bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns); 103 103 private: 104 104 bool IsSymmetric; -
src/periodentafel.cpp
r51c910 r36ec71 92 92 * \return pointer to element 93 93 */ 94 element * periodentafel::FindElement(c har *shorthand) const94 element * periodentafel::FindElement(const char *shorthand) const 95 95 { 96 96 element *walker = periodentafel::start; … … 162 162 * \param *path to to standard file names 163 163 */ 164 bool periodentafel::LoadPeriodentafel(c har *path)164 bool periodentafel::LoadPeriodentafel(const char *path) 165 165 { 166 166 ifstream infile; … … 301 301 /** Stores element list to file. 302 302 */ 303 bool periodentafel::StorePeriodentafel(c har *path) const303 bool periodentafel::StorePeriodentafel(const char *path) const 304 304 { 305 305 bool result = true; -
src/periodentafel.hpp
r51c910 r36ec71 21 21 */ 22 22 class element { 23 24 double mass;//!< mass in g/mol25 double CovalentRadius;//!< covalent radius26 double VanDerWaalsRadius;//!< can-der-Waals radius27 int Z;//!< atomic number28 char name[64];//!< atom name, i.e. "Hydrogren"29 30 char period[8];//!< period: n quantum number31 char group[8];//!< group: l quantum number32 char block[8];//!< block: l quantum number33 element *previous;//!< previous item in list34 element *next;//!< next element in list35 int *sort;//!< sorc criteria36 int No;//!< number of element set on periodentafel::Output()37 double Valence;//!< number of valence electrons for this element38 int NoValenceOrbitals;//!< number of valence orbitals, used for determining bond degree in molecule::CreateConnectmatrix()39 double HBondDistance[NDIM]; //!< distance in Angstrom of this element to hydrogen(for single, double and triple bonds)40 double HBondAngle[NDIM];//!< typical angle for one, two, three bonded hydrogen (in degrees)23 public: 24 double mass; //!< mass in g/mol 25 double CovalentRadius; //!< covalent radius 26 double VanDerWaalsRadius; //!< can-der-Waals radius 27 int Z; //!< atomic number 28 char name[64]; //!< atom name, i.e. "Hydrogren" 29 char symbol[3]; //!< short form of the atom, i.e. "H" 30 char period[8]; //!< period: n quantum number 31 char group[8]; //!< group: l quantum number 32 char block[8]; //!< block: l quantum number 33 element *previous; //!< previous item in list 34 element *next; //!< next element in list 35 int *sort; //!< sorc criteria 36 int No; //!< number of element set on periodentafel::Output() 37 double Valence; //!< number of valence electrons for this element 38 int NoValenceOrbitals; //!< number of valence orbitals, used for determining bond degree in molecule::CreateConnectmatrix() 39 double HBondDistance[NDIM]; //!< distance in Angstrom of this element to hydrogen (for single, double and triple bonds) 40 double HBondAngle[NDIM]; //!< typical angle for one, two, three bonded hydrogen (in degrees) 41 41 42 43 42 element(); 43 ~element(); 44 44 45 46 47 45 //> print element entries to screen 46 bool Output(ofstream *out) const; 47 bool Checkout(ofstream *out, const int No, const int NoOfAtoms) const; 48 48 49 49 private: 50 50 }; 51 51 … … 53 53 */ 54 54 class periodentafel { 55 56 57 element *end;//!< end of element list58 59 55 public: 56 element *start; //!< start of element list 57 element *end; //!< end of element list 58 char header1[MAXSTRINGSIZE]; //!< store first header line 59 char header2[MAXSTRINGSIZE]; //!< store second header line 60 60 61 62 61 periodentafel(); 62 ~periodentafel(); 63 63 64 65 66 67 68 element * FindElement(char *shorthand) const;69 70 71 72 bool LoadPeriodentafel(char *path = NULL);73 bool StorePeriodentafel(char *path = NULL) const;64 bool AddElement(element *pointer); 65 bool RemoveElement(element *pointer); 66 bool CleanupPeriodtable(); 67 element * FindElement(int Z); 68 element * FindElement(const char *shorthand) const; 69 element * AskElement(); 70 bool Output(ofstream *output) const; 71 bool Checkout(ofstream *output, const int *checkliste) const; 72 bool LoadPeriodentafel(const char *path); 73 bool StorePeriodentafel(const char *path) const; 74 74 75 75 private: 76 76 }; 77 77 -
src/stackclass.hpp
r51c910 r36ec71 10 10 */ 11 11 template <typename T> class StackClass { 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 T *StackList;//!< the list containing the atom pointers29 int EntryCount;//!< number of entries in the stack30 int CurrentLastEntry;//!< Current last entry (newest item on stack)31 int CurrentFirstEntry;//!< Current first entry (oldest item on stack)32 int NextFreeField;//!< Current index of next free field12 public: 13 StackClass<T>(int dimension); 14 ~StackClass<T>(); 15 16 bool Push(T object); 17 T PopFirst(); 18 T PopLast(); 19 bool RemoveItem(T ptr); 20 void ClearStack(); 21 bool IsEmpty(); 22 bool IsFull(); 23 int ItemCount(); 24 void Output(ofstream *out) const; 25 void TestImplementation(ofstream *out, T test); 26 27 private: 28 T *StackList; //!< the list containing the atom pointers 29 int EntryCount; //!< number of entries in the stack 30 int CurrentLastEntry; //!< Current last entry (newest item on stack) 31 int CurrentFirstEntry; //!< Current first entry (oldest item on stack) 32 int NextFreeField; //!< Current index of next free field 33 33 }; 34 34 … … 37 37 template <typename T> StackClass<T>::StackClass(int dimension) 38 38 { 39 40 41 42 43 39 CurrentLastEntry = 0; 40 CurrentFirstEntry = 0; 41 NextFreeField = 0; 42 EntryCount = dimension; 43 StackList = (T *) Malloc(sizeof(T)*EntryCount, "StackClass::StackClass: **StackList"); 44 44 }; 45 45 … … 48 48 template <typename T> StackClass<T>::~StackClass() 49 49 { 50 50 Free((void **)&StackList, "StackClass::StackClass: **StackList"); 51 51 }; 52 52 … … 57 57 template <typename T> bool StackClass<T>::Push(T object) 58 58 { 59 if (!IsFull()) {// check whether free field is really not occupied60 61 62 63 64 65 66 67 59 if (!IsFull()) { // check whether free field is really not occupied 60 StackList[NextFreeField] = object; 61 CurrentLastEntry = NextFreeField; 62 NextFreeField = (NextFreeField + 1) % EntryCount; // step on to next free field 63 return true; 64 } else { 65 cerr << "ERROR: Stack is full, " << "Stack: CurrentLastEntry " << CurrentLastEntry<< "\tCurrentFirstEntry " << CurrentFirstEntry << "\tNextFreeField " << NextFreeField << "\tEntryCount " << EntryCount << "!" << endl; 66 return false; 67 } 68 68 }; 69 69 … … 74 74 template <typename T> T StackClass<T>::PopFirst() 75 75 { 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 76 T Walker = NULL; 77 if (!IsEmpty()) { 78 Walker = StackList[CurrentFirstEntry]; 79 if (Walker == NULL) 80 cerr << "ERROR: Stack's field is empty!" << endl; 81 StackList[CurrentFirstEntry] = NULL; 82 if (CurrentFirstEntry != CurrentLastEntry) { // hasn't last item been popped as well? 83 CurrentFirstEntry = (CurrentFirstEntry + 1) % EntryCount; // step back from current free field to last used (somehow modulo does not work in -1) 84 } else { 85 CurrentFirstEntry = (CurrentFirstEntry + 1) % EntryCount; // step back from current free field to last used (somehow modulo does not work in -1) 86 CurrentLastEntry = CurrentFirstEntry; 87 } 88 } else 89 cerr << "ERROR: Stack is empty!" << endl; 90 return Walker; 91 91 }; 92 92 … … 97 97 template <typename T> T StackClass<T>::PopLast() 98 98 { 99 100 101 102 103 104 105 106 if (CurrentLastEntry != CurrentFirstEntry)// has there been more than one item on stack107 108 109 110 111 99 T Walker = NULL; 100 if (!IsEmpty()) { 101 Walker = StackList[CurrentLastEntry]; 102 StackList[CurrentLastEntry] = NULL; 103 if (Walker == NULL) 104 cerr << "ERROR: Stack's field is empty!" << endl; 105 NextFreeField = CurrentLastEntry; 106 if (CurrentLastEntry != CurrentFirstEntry) // has there been more than one item on stack 107 CurrentLastEntry = (CurrentLastEntry + (EntryCount-1)) % EntryCount; // step back from current free field to last (modulo does not work in -1, thus go EntryCount-1 instead) 108 } else { 109 cerr << "ERROR: Stack is empty!" << endl; 110 } 111 return Walker; 112 112 }; 113 113 … … 120 120 template <typename T> bool StackClass<T>::RemoveItem(T ptr) 121 121 { 122 123 124 125 126 127 if (StackList[i] == ptr) {// if item found, remove128 129 130 131 132 if ((found) && (StackList[i] != NULL)) {// means we have to shift (and not the removed item)133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 if (CurrentLastEntry != CurrentFirstEntry)// has there been more than one item on stack148 149 150 122 bool found = false; 123 cout << Verbose(5) << "First " << CurrentFirstEntry<< "\tLast " << CurrentLastEntry<< "\tNext " << NextFreeField<< "\tCount " << EntryCount<< "." << endl; 124 int i=CurrentFirstEntry; 125 if (!IsEmpty()) 126 do { 127 if (StackList[i] == ptr) { // if item found, remove 128 cout << Verbose(5) << "Item " << *ptr << " was number " << i << " on stack, removing it." << endl; 129 found = true; 130 StackList[i] = NULL; 131 } 132 if ((found) && (StackList[i] != NULL)) { // means we have to shift (and not the removed item) 133 if (i == 0) { // we are down to first item in stack, have to put onto last item 134 cout << Verbose(5) << "Shifting item 0 to place " << EntryCount-1 << "." << endl; 135 StackList[EntryCount-1] = StackList[0]; 136 } else { 137 cout << Verbose(5) << "Shifting item " << i << " to place " << i-1 << "." << endl; 138 StackList[i-1] = StackList[i]; 139 } 140 } 141 i=((i + 1) % EntryCount); // step on 142 } while (i!=NextFreeField); 143 else 144 cerr << "ERROR: Stack is already empty!" << endl; 145 if (found) { 146 NextFreeField = CurrentLastEntry; 147 if (CurrentLastEntry != CurrentFirstEntry) // has there been more than one item on stack 148 CurrentLastEntry = (CurrentLastEntry + (EntryCount-1)) % EntryCount; 149 } 150 return found; 151 151 }; 152 152 … … 158 158 template <typename T> void StackClass<T>::TestImplementation(ofstream *out, T test) 159 159 { 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 *out << "Removing first element..." << endl;180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 160 T Walker = test; 161 *out << Verbose(1) << "Testing the snake stack..." << endl; 162 for (int i=0;i<EntryCount;i++) { 163 *out << Verbose(2) << "Filling " << i << "th element of stack." << endl; 164 Push(Walker); 165 Walker=Walker->next; 166 } 167 *out << endl; 168 Output(out); 169 if (IsFull()) 170 *out << "Stack is full as supposed to be!" << endl; 171 else 172 *out << "ERROR: Stack is not as full as supposed to be!" << endl; 173 //if (StackList[(EntryCount+1)/2] != NULL) { 174 *out << "Removing element in the middle ..." << endl; 175 RemoveItem(StackList[(EntryCount+1)/2]); 176 Output(out); 177 //} 178 //if (StackList[CurrentFirstEntry] != NULL) { 179 *out << "Removing first element ..." << endl; 180 RemoveItem(StackList[CurrentFirstEntry]); 181 Output(out); 182 //} 183 //if (StackList[CurrentLastEntry] != NULL) { 184 *out << "Removing last element ..." << endl; 185 RemoveItem(StackList[CurrentLastEntry]); 186 Output(out); 187 //} 188 *out << "Clearing stack ... " << endl; 189 ClearStack(); 190 Output(out); 191 if (IsEmpty()) 192 *out << "Stack is empty as supposed to be!" << endl; 193 else 194 *out << "ERROR: Stack is not as empty as supposed to be!" << endl; 195 *out << "done." << endl; 196 196 }; 197 197 … … 201 201 template <typename T> void StackClass<T>::Output(ofstream *out) const 202 202 { 203 204 205 206 207 208 if(i == CurrentLastEntry)209 210 if (i ==NextFreeField)211 212 213 214 203 *out << "Contents of Stack: "; 204 for(int i=0;i<EntryCount;i++) { 205 *out << "\t"; 206 if (i == CurrentFirstEntry) 207 *out << " 1"; 208 if (i == CurrentLastEntry) 209 *out << " "<< EntryCount; 210 if (i == NextFreeField) 211 *out << " F"; 212 *out << ": " << StackList[i]; 213 } 214 *out << endl; 215 215 }; 216 216 … … 222 222 template <typename T> bool StackClass<T>::IsEmpty() 223 223 { 224 224 return((NextFreeField == CurrentFirstEntry) && (CurrentLastEntry == CurrentFirstEntry)); 225 225 }; 226 226 … … 232 232 template <typename T> bool StackClass<T>::IsFull() 233 233 { 234 234 return((NextFreeField == CurrentFirstEntry) && (CurrentLastEntry != CurrentFirstEntry)); 235 235 }; 236 236 … … 242 242 template <typename T> int StackClass<T>::ItemCount() 243 243 { 244 245 244 //cout << "Stack: CurrentLastEntry " << CurrentLastEntry<< "\tCurrentFirstEntry " << CurrentFirstEntry << "\tEntryCount " << EntryCount << "." << endl; 245 return((NextFreeField + (EntryCount - CurrentFirstEntry)) % EntryCount); 246 246 }; 247 247 … … 251 251 template <typename T> void StackClass<T>::ClearStack() 252 252 { 253 254 255 256 257 253 for(int i=EntryCount; i--;) 254 StackList[i] = NULL; 255 CurrentFirstEntry = 0; 256 CurrentLastEntry = 0; 257 NextFreeField = 0; 258 258 }; 259 259 -
src/vector.cpp
r51c910 r36ec71 28 28 double Vector::DistanceSquared(const Vector *y) const 29 29 { 30 31 32 33 30 double res = 0.; 31 for (int i=NDIM;i--;) 32 res += (x[i]-y->x[i])*(x[i]-y->x[i]); 33 return (res); 34 34 }; 35 35 … … 40 40 double Vector::Distance(const Vector *y) const 41 41 { 42 43 44 45 42 double res = 0.; 43 for (int i=NDIM;i--;) 44 res += (x[i]-y->x[i])*(x[i]-y->x[i]); 45 return (sqrt(res)); 46 46 }; 47 47 … … 53 53 double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const 54 54 { 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 55 double res = Distance(y), tmp, matrix[NDIM*NDIM]; 56 Vector Shiftedy, TranslationVector; 57 int N[NDIM]; 58 matrix[0] = cell_size[0]; 59 matrix[1] = cell_size[1]; 60 matrix[2] = cell_size[3]; 61 matrix[3] = cell_size[1]; 62 matrix[4] = cell_size[2]; 63 matrix[5] = cell_size[4]; 64 matrix[6] = cell_size[3]; 65 matrix[7] = cell_size[4]; 66 matrix[8] = cell_size[5]; 67 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells 68 for (N[0]=-1;N[0]<=1;N[0]++) 69 for (N[1]=-1;N[1]<=1;N[1]++) 70 for (N[2]=-1;N[2]<=1;N[2]++) { 71 // create the translation vector 72 TranslationVector.Zero(); 73 for (int i=NDIM;i--;) 74 TranslationVector.x[i] = (double)N[i]; 75 TranslationVector.MatrixMultiplication(matrix); 76 // add onto the original vector to compare with 77 Shiftedy.CopyVector(y); 78 Shiftedy.AddVector(&TranslationVector); 79 // get distance and compare with minimum so far 80 tmp = Distance(&Shiftedy); 81 if (tmp < res) res = tmp; 82 } 83 return (res); 84 84 }; 85 85 … … 91 91 double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const 92 92 { 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 93 double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM]; 94 Vector Shiftedy, TranslationVector; 95 int N[NDIM]; 96 matrix[0] = cell_size[0]; 97 matrix[1] = cell_size[1]; 98 matrix[2] = cell_size[3]; 99 matrix[3] = cell_size[1]; 100 matrix[4] = cell_size[2]; 101 matrix[5] = cell_size[4]; 102 matrix[6] = cell_size[3]; 103 matrix[7] = cell_size[4]; 104 matrix[8] = cell_size[5]; 105 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells 106 for (N[0]=-1;N[0]<=1;N[0]++) 107 for (N[1]=-1;N[1]<=1;N[1]++) 108 for (N[2]=-1;N[2]<=1;N[2]++) { 109 // create the translation vector 110 TranslationVector.Zero(); 111 for (int i=NDIM;i--;) 112 TranslationVector.x[i] = (double)N[i]; 113 TranslationVector.MatrixMultiplication(matrix); 114 // add onto the original vector to compare with 115 Shiftedy.CopyVector(y); 116 Shiftedy.AddVector(&TranslationVector); 117 // get distance and compare with minimum so far 118 tmp = DistanceSquared(&Shiftedy); 119 if (tmp < res) res = tmp; 120 } 121 return (res); 122 122 }; 123 123 … … 128 128 void Vector::KeepPeriodic(ofstream *out, double *matrix) 129 129 { 130 // 131 // 132 133 134 // 135 // 136 // 137 // 138 139 140 141 if (TestVector.x[i] < 0) {// get every coefficient into the interval [0,1)142 143 144 145 146 147 148 149 // 150 // 151 // 152 // 130 // int N[NDIM]; 131 // bool flag = false; 132 //vector Shifted, TranslationVector; 133 Vector TestVector; 134 // *out << Verbose(1) << "Begin of KeepPeriodic." << endl; 135 // *out << Verbose(2) << "Vector is: "; 136 // Output(out); 137 // *out << endl; 138 TestVector.CopyVector(this); 139 TestVector.InverseMatrixMultiplication(matrix); 140 for(int i=NDIM;i--;) { // correct periodically 141 if (TestVector.x[i] < 0) { // get every coefficient into the interval [0,1) 142 TestVector.x[i] += ceil(TestVector.x[i]); 143 } else { 144 TestVector.x[i] -= floor(TestVector.x[i]); 145 } 146 } 147 TestVector.MatrixMultiplication(matrix); 148 CopyVector(&TestVector); 149 // *out << Verbose(2) << "New corrected vector is: "; 150 // Output(out); 151 // *out << endl; 152 // *out << Verbose(1) << "End of KeepPeriodic." << endl; 153 153 }; 154 154 … … 159 159 double Vector::ScalarProduct(const Vector *y) const 160 160 { 161 162 163 164 161 double res = 0.; 162 for (int i=NDIM;i--;) 163 res += x[i]*y->x[i]; 164 return (res); 165 165 }; 166 166 167 167 168 168 /** Calculates VectorProduct between this and another vector. 169 * 170 * 171 * 172 * 169 * -# returns the Product in place of vector from which it was initiated 170 * -# ATTENTION: Only three dim. 171 * \param *y array to vector with which to calculate crossproduct 172 * \return \f$ x \times y \f& 173 173 */ 174 174 void Vector::VectorProduct(const Vector *y) 175 175 { 176 177 178 179 180 176 Vector tmp; 177 tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]); 178 tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]); 179 tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]); 180 this->CopyVector(&tmp); 181 181 182 182 }; … … 189 189 void Vector::ProjectOntoPlane(const Vector *y) 190 190 { 191 192 193 194 195 191 Vector tmp; 192 tmp.CopyVector(y); 193 tmp.Normalize(); 194 tmp.Scale(ScalarProduct(&tmp)); 195 this->SubtractVector(&tmp); 196 196 }; 197 197 … … 202 202 double Vector::Projection(const Vector *y) const 203 203 { 204 204 return (ScalarProduct(y)); 205 205 }; 206 206 … … 210 210 double Vector::Norm() const 211 211 { 212 213 214 215 212 double res = 0.; 213 for (int i=NDIM;i--;) 214 res += this->x[i]*this->x[i]; 215 return (sqrt(res)); 216 216 }; 217 217 … … 220 220 void Vector::Normalize() 221 221 { 222 223 224 225 226 227 222 double res = 0.; 223 for (int i=NDIM;i--;) 224 res += this->x[i]*this->x[i]; 225 if (fabs(res) > MYEPSILON) 226 res = 1./sqrt(res); 227 Scale(&res); 228 228 }; 229 229 … … 232 232 void Vector::Zero() 233 233 { 234 235 234 for (int i=NDIM;i--;) 235 this->x[i] = 0.; 236 236 }; 237 237 … … 240 240 void Vector::One(double one) 241 241 { 242 243 242 for (int i=NDIM;i--;) 243 this->x[i] = one; 244 244 }; 245 245 … … 248 248 void Vector::Init(double x1, double x2, double x3) 249 249 { 250 251 252 250 x[0] = x1; 251 x[1] = x2; 252 x[2] = x3; 253 253 }; 254 254 … … 259 259 double Vector::Angle(const Vector *y) const 260 260 { 261 return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 261 double angle = this->ScalarProduct(y)/Norm()/y->Norm(); 262 // -1-MYEPSILON occured due to numerical imprecision, catch ... 263 //cout << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl; 264 if (angle < -1) 265 angle = -1; 266 if (angle > 1) 267 angle = 1; 268 return acos(angle); 262 269 }; 263 270 … … 268 275 void Vector::RotateVector(const Vector *axis, const double alpha) 269 276 { 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 277 Vector a,y; 278 // normalise this vector with respect to axis 279 a.CopyVector(this); 280 a.Scale(Projection(axis)); 281 SubtractVector(&a); 282 // construct normal vector 283 y.MakeNormalVector(axis,this); 284 y.Scale(Norm()); 285 // scale normal vector by sine and this vector by cosine 286 y.Scale(sin(alpha)); 287 Scale(cos(alpha)); 288 // add scaled normal vector onto this vector 289 AddVector(&y); 290 // add part in axis direction 291 AddVector(&a); 285 292 }; 286 293 … … 292 299 Vector& operator+=(Vector& a, const Vector& b) 293 300 { 294 295 301 a.AddVector(&b); 302 return a; 296 303 }; 297 304 /** factor each component of \a a times a double \a m. … … 302 309 Vector& operator*=(Vector& a, const double m) 303 310 { 304 305 306 }; 307 308 /** Sums two vectors \a 311 a.Scale(m); 312 return a; 313 }; 314 315 /** Sums two vectors \a and \b component-wise. 309 316 * \param a first vector 310 317 * \param b second vector … … 313 320 Vector& operator+(const Vector& a, const Vector& b) 314 321 { 315 316 317 318 322 Vector *x = new Vector; 323 x->CopyVector(&a); 324 x->AddVector(&b); 325 return *x; 319 326 }; 320 327 … … 326 333 Vector& operator*(const Vector& a, const double m) 327 334 { 328 329 330 331 335 Vector *x = new Vector; 336 x->CopyVector(&a); 337 x->Scale(m); 338 return *x; 332 339 }; 333 340 … … 338 345 bool Vector::Output(ofstream *out) const 339 346 { 340 341 342 343 344 345 346 347 348 349 350 347 if (out != NULL) { 348 *out << "("; 349 for (int i=0;i<NDIM;i++) { 350 *out << x[i]; 351 if (i != 2) 352 *out << ","; 353 } 354 *out << ")"; 355 return true; 356 } else 357 return false; 351 358 }; 352 359 353 360 ostream& operator<<(ostream& ost,Vector& m) 354 361 { 355 356 357 358 359 360 361 362 362 ost << "("; 363 for (int i=0;i<NDIM;i++) { 364 ost << m.x[i]; 365 if (i != 2) 366 ost << ","; 367 } 368 ost << ")"; 369 return ost; 363 370 }; 364 371 … … 368 375 void Vector::Scale(double **factor) 369 376 { 370 371 377 for (int i=NDIM;i--;) 378 x[i] *= (*factor)[i]; 372 379 }; 373 380 374 381 void Vector::Scale(double *factor) 375 382 { 376 377 383 for (int i=NDIM;i--;) 384 x[i] *= *factor; 378 385 }; 379 386 380 387 void Vector::Scale(double factor) 381 388 { 382 383 389 for (int i=NDIM;i--;) 390 x[i] *= factor; 384 391 }; 385 392 … … 389 396 void Vector::Translate(const Vector *trans) 390 397 { 391 392 398 for (int i=NDIM;i--;) 399 x[i] += trans->x[i]; 393 400 }; 394 401 … … 398 405 void Vector::MatrixMultiplication(double *M) 399 406 { 400 Vector C; 401 // do the matrix multiplication 402 C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2]; 403 C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2]; 404 C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2]; 405 // transfer the result into this 406 for (int i=NDIM;i--;) 407 x[i] = C.x[i]; 407 Vector C; 408 // do the matrix multiplication 409 C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2]; 410 C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2]; 411 C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2]; 412 // transfer the result into this 413 for (int i=NDIM;i--;) 414 x[i] = C.x[i]; 415 }; 416 417 /** Calculate the inverse of a 3x3 matrix. 418 * \param *matrix NDIM_NDIM array 419 */ 420 double * Vector::InverseMatrix(double *A) 421 { 422 double *B = (double *) Malloc(sizeof(double)*NDIM*NDIM, "Vector::InverseMatrix: *B"); 423 double detA = RDET3(A); 424 double detAReci; 425 426 for (int i=0;i<NDIM*NDIM;++i) 427 B[i] = 0.; 428 // calculate the inverse B 429 if (fabs(detA) > MYEPSILON) {; // RDET3(A) yields precisely zero if A irregular 430 detAReci = 1./detA; 431 B[0] = detAReci*RDET2(A[4],A[5],A[7],A[8]); // A_11 432 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]); // A_12 433 B[2] = detAReci*RDET2(A[1],A[2],A[4],A[5]); // A_13 434 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]); // A_21 435 B[4] = detAReci*RDET2(A[0],A[2],A[6],A[8]); // A_22 436 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]); // A_23 437 B[6] = detAReci*RDET2(A[3],A[4],A[6],A[7]); // A_31 438 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]); // A_32 439 B[8] = detAReci*RDET2(A[0],A[1],A[3],A[4]); // A_33 440 } 441 return B; 408 442 }; 409 443 … … 413 447 void Vector::InverseMatrixMultiplication(double *A) 414 448 { 415 416 417 418 419 420 421 if (fabs(detA) > MYEPSILON) {;// RDET3(A) yields precisely zero if A irregular422 423 B[0] = detAReci*RDET2(A[4],A[5],A[7],A[8]);// A_11424 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);// A_12425 B[2] = detAReci*RDET2(A[1],A[2],A[4],A[5]);// A_13426 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);// A_21427 B[4] = detAReci*RDET2(A[0],A[2],A[6],A[8]);// A_22428 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);// A_23429 B[6] = detAReci*RDET2(A[3],A[4],A[6],A[7]);// A_31430 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);// A_32431 B[8] = detAReci*RDET2(A[0],A[1],A[3],A[4]);// A_33432 433 434 435 436 437 438 439 440 441 442 449 Vector C; 450 double B[NDIM*NDIM]; 451 double detA = RDET3(A); 452 double detAReci; 453 454 // calculate the inverse B 455 if (fabs(detA) > MYEPSILON) {; // RDET3(A) yields precisely zero if A irregular 456 detAReci = 1./detA; 457 B[0] = detAReci*RDET2(A[4],A[5],A[7],A[8]); // A_11 458 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]); // A_12 459 B[2] = detAReci*RDET2(A[1],A[2],A[4],A[5]); // A_13 460 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]); // A_21 461 B[4] = detAReci*RDET2(A[0],A[2],A[6],A[8]); // A_22 462 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]); // A_23 463 B[6] = detAReci*RDET2(A[3],A[4],A[6],A[7]); // A_31 464 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]); // A_32 465 B[8] = detAReci*RDET2(A[0],A[1],A[3],A[4]); // A_33 466 467 // do the matrix multiplication 468 C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2]; 469 C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2]; 470 C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2]; 471 // transfer the result into this 472 for (int i=NDIM;i--;) 473 x[i] = C.x[i]; 474 } else { 475 cerr << "ERROR: inverse of matrix does not exists!" << endl; 476 } 443 477 }; 444 478 … … 453 487 void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors) 454 488 { 455 456 489 for(int i=NDIM;i--;) 490 x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i]; 457 491 }; 458 492 … … 462 496 void Vector::Mirror(const Vector *n) 463 497 { 464 465 projection = ScalarProduct(n)/n->ScalarProduct(n);// remove constancy from n (keep as logical one)466 467 468 469 470 471 472 473 474 498 double projection; 499 projection = ScalarProduct(n)/n->ScalarProduct(n); // remove constancy from n (keep as logical one) 500 // withdraw projected vector twice from original one 501 cout << Verbose(1) << "Vector: "; 502 Output((ofstream *)&cout); 503 cout << "\t"; 504 for (int i=NDIM;i--;) 505 x[i] -= 2.*projection*n->x[i]; 506 cout << "Projected vector: "; 507 Output((ofstream *)&cout); 508 cout << endl; 475 509 }; 476 510 … … 484 518 bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3) 485 519 { 486 487 488 489 490 491 492 493 494 495 496 // 497 // 498 // 499 // 500 // 501 // 502 503 504 505 506 507 508 520 Vector x1, x2; 521 522 x1.CopyVector(y1); 523 x1.SubtractVector(y2); 524 x2.CopyVector(y3); 525 x2.SubtractVector(y2); 526 if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) { 527 cout << Verbose(4) << "Given vectors are linear dependent." << endl; 528 return false; 529 } 530 // cout << Verbose(4) << "relative, first plane coordinates:"; 531 // x1.Output((ofstream *)&cout); 532 // cout << endl; 533 // cout << Verbose(4) << "second plane coordinates:"; 534 // x2.Output((ofstream *)&cout); 535 // cout << endl; 536 537 this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]); 538 this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]); 539 this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]); 540 Normalize(); 541 542 return true; 509 543 }; 510 544 … … 520 554 bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2) 521 555 { 522 523 524 525 526 527 528 529 530 // 531 // 532 // 533 // 534 // 535 // 536 537 538 539 540 541 542 556 Vector x1,x2; 557 x1.CopyVector(y1); 558 x2.CopyVector(y2); 559 Zero(); 560 if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) { 561 cout << Verbose(4) << "Given vectors are linear dependent." << endl; 562 return false; 563 } 564 // cout << Verbose(4) << "relative, first plane coordinates:"; 565 // x1.Output((ofstream *)&cout); 566 // cout << endl; 567 // cout << Verbose(4) << "second plane coordinates:"; 568 // x2.Output((ofstream *)&cout); 569 // cout << endl; 570 571 this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]); 572 this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]); 573 this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]); 574 Normalize(); 575 576 return true; 543 577 }; 544 578 … … 550 584 bool Vector::MakeNormalVector(const Vector *y1) 551 585 { 552 553 554 555 556 557 558 559 560 586 bool result = false; 587 Vector x1; 588 x1.CopyVector(y1); 589 x1.Scale(x1.Projection(this)); 590 SubtractVector(&x1); 591 for (int i=NDIM;i--;) 592 result = result || (fabs(x[i]) > MYEPSILON); 593 594 return result; 561 595 }; 562 596 … … 569 603 bool Vector::GetOneNormalVector(const Vector *GivenVector) 570 604 { 571 572 int Last = 0;// count the number of non-zero entries in vector573 int j;// loop variables574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 case 3:// threecomponent system589 case 2:// two component system590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 605 int Components[NDIM]; // contains indices of non-zero components 606 int Last = 0; // count the number of non-zero entries in vector 607 int j; // loop variables 608 double norm; 609 610 cout << Verbose(4); 611 GivenVector->Output((ofstream *)&cout); 612 cout << endl; 613 for (j=NDIM;j--;) 614 Components[j] = -1; 615 // find two components != 0 616 for (j=0;j<NDIM;j++) 617 if (fabs(GivenVector->x[j]) > MYEPSILON) 618 Components[Last++] = j; 619 cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl; 620 621 switch(Last) { 622 case 3: // threecomponent system 623 case 2: // two component system 624 norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]])); 625 x[Components[2]] = 0.; 626 // in skp both remaining parts shall become zero but with opposite sign and third is zero 627 x[Components[1]] = -1./GivenVector->x[Components[1]] / norm; 628 x[Components[0]] = 1./GivenVector->x[Components[0]] / norm; 629 return true; 630 break; 631 case 1: // one component system 632 // set sole non-zero component to 0, and one of the other zero component pendants to 1 633 x[(Components[0]+2)%NDIM] = 0.; 634 x[(Components[0]+1)%NDIM] = 1.; 635 x[Components[0]] = 0.; 636 return true; 637 break; 638 default: 639 return false; 640 } 607 641 }; 608 642 … … 615 649 double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C) 616 650 { 617 // 618 // 619 // 620 // 621 // 622 651 // cout << Verbose(3) << "For comparison: "; 652 // cout << "A " << A->Projection(this) << "\t"; 653 // cout << "B " << B->Projection(this) << "\t"; 654 // cout << "C " << C->Projection(this) << "\t"; 655 // cout << endl; 656 return A->Projection(this); 623 657 }; 624 658 … … 630 664 bool Vector::LSQdistance(Vector **vectors, int num) 631 665 { 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 666 int j; 667 668 for (j=0;j<num;j++) { 669 cout << Verbose(1) << j << "th atom's vector: "; 670 (vectors[j])->Output((ofstream *)&cout); 671 cout << endl; 672 } 673 674 int np = 3; 675 struct LSQ_params par; 676 677 const gsl_multimin_fminimizer_type *T = 678 gsl_multimin_fminimizer_nmsimplex; 679 gsl_multimin_fminimizer *s = NULL; 680 gsl_vector *ss, *y; 681 gsl_multimin_function minex_func; 682 683 size_t iter = 0, i; 684 int status; 685 double size; 686 687 /* Initial vertex size vector */ 688 ss = gsl_vector_alloc (np); 689 y = gsl_vector_alloc (np); 690 691 /* Set all step sizes to 1 */ 692 gsl_vector_set_all (ss, 1.0); 693 694 /* Starting point */ 695 par.vectors = vectors; 696 par.num = num; 697 698 for (i=NDIM;i--;) 699 gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.); 700 701 /* Initialize method and iterate */ 702 minex_func.f = &LSQ; 703 minex_func.n = np; 704 minex_func.params = (void *)∥ 705 706 s = gsl_multimin_fminimizer_alloc (T, np); 707 gsl_multimin_fminimizer_set (s, &minex_func, y, ss); 708 709 do 710 { 711 iter++; 712 status = gsl_multimin_fminimizer_iterate(s); 713 714 if (status) 715 break; 716 717 size = gsl_multimin_fminimizer_size (s); 718 status = gsl_multimin_test_size (size, 1e-2); 719 720 if (status == GSL_SUCCESS) 721 { 722 printf ("converged to minimum at\n"); 723 } 724 725 printf ("%5d ", (int)iter); 726 for (i = 0; i < (size_t)np; i++) 727 { 728 printf ("%10.3e ", gsl_vector_get (s->x, i)); 729 } 730 printf ("f() = %7.3f size = %.3f\n", s->fval, size); 731 } 732 while (status == GSL_CONTINUE && iter < 100); 733 734 for (i=(size_t)np;i--;) 735 this->x[i] = gsl_vector_get(s->x, i); 736 gsl_vector_free(y); 737 gsl_vector_free(ss); 738 gsl_multimin_fminimizer_free (s); 739 740 return true; 707 741 }; 708 742 … … 712 746 void Vector::AddVector(const Vector *y) 713 747 { 714 715 748 for (int i=NDIM;i--;) 749 this->x[i] += y->x[i]; 716 750 } 717 751 … … 721 755 void Vector::SubtractVector(const Vector *y) 722 756 { 723 724 757 for (int i=NDIM;i--;) 758 this->x[i] -= y->x[i]; 725 759 } 726 760 … … 730 764 void Vector::CopyVector(const Vector *y) 731 765 { 732 733 766 for (int i=NDIM;i--;) 767 this->x[i] = y->x[i]; 734 768 } 735 769 … … 741 775 void Vector::AskPosition(double *cell_size, bool check) 742 776 { 743 744 745 746 747 748 749 750 751 777 char coords[3] = {'x','y','z'}; 778 int j = -1; 779 for (int i=0;i<3;i++) { 780 j += i+1; 781 do { 782 cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: "; 783 cin >> x[i]; 784 } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check)); 785 } 752 786 }; 753 787 … … 769 803 bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c) 770 804 { 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 x[0] =A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 }; 805 double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C; 806 double ang; // angle on testing 807 double sign[3]; 808 int i,j,k; 809 A = cos(alpha) * x1->Norm() * c; 810 B1 = cos(beta + M_PI/2.) * y->Norm() * c; 811 B2 = cos(beta) * x2->Norm() * c; 812 C = c * c; 813 cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl; 814 int flag = 0; 815 if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping 816 if (fabs(x1->x[1]) > MYEPSILON) { 817 flag = 1; 818 } else if (fabs(x1->x[2]) > MYEPSILON) { 819 flag = 2; 820 } else { 821 return false; 822 } 823 } 824 switch (flag) { 825 default: 826 case 0: 827 break; 828 case 2: 829 flip(&x1->x[0],&x1->x[1]); 830 flip(&x2->x[0],&x2->x[1]); 831 flip(&y->x[0],&y->x[1]); 832 //flip(&x[0],&x[1]); 833 flip(&x1->x[1],&x1->x[2]); 834 flip(&x2->x[1],&x2->x[2]); 835 flip(&y->x[1],&y->x[2]); 836 //flip(&x[1],&x[2]); 837 case 1: 838 flip(&x1->x[0],&x1->x[1]); 839 flip(&x2->x[0],&x2->x[1]); 840 flip(&y->x[0],&y->x[1]); 841 //flip(&x[0],&x[1]); 842 flip(&x1->x[1],&x1->x[2]); 843 flip(&x2->x[1],&x2->x[2]); 844 flip(&y->x[1],&y->x[2]); 845 //flip(&x[1],&x[2]); 846 break; 847 } 848 // now comes the case system 849 D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1]; 850 D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 851 D3 = y->x[0]/x1->x[0]*A-B1; 852 cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n"; 853 if (fabs(D1) < MYEPSILON) { 854 cout << Verbose(2) << "D1 == 0!\n"; 855 if (fabs(D2) > MYEPSILON) { 856 cout << Verbose(3) << "D2 != 0!\n"; 857 x[2] = -D3/D2; 858 E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2; 859 E2 = -x1->x[1]/x1->x[0]; 860 cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n"; 861 F1 = E1*E1 + 1.; 862 F2 = -E1*E2; 863 F3 = E1*E1 + D3*D3/(D2*D2) - C; 864 cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n"; 865 if (fabs(F1) < MYEPSILON) { 866 cout << Verbose(4) << "F1 == 0!\n"; 867 cout << Verbose(4) << "Gleichungssystem linear\n"; 868 x[1] = F3/(2.*F2); 869 } else { 870 p = F2/F1; 871 q = p*p - F3/F1; 872 cout << Verbose(4) << "p " << p << "\tq " << q << endl; 873 if (q < 0) { 874 cout << Verbose(4) << "q < 0" << endl; 875 return false; 876 } 877 x[1] = p + sqrt(q); 878 } 879 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2]; 880 } else { 881 cout << Verbose(2) << "Gleichungssystem unterbestimmt\n"; 882 return false; 883 } 884 } else { 885 E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1; 886 E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2]; 887 cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n"; 888 F1 = E2*E2 + D2*D2/(D1*D1) + 1.; 889 F2 = -(E1*E2 + D2*D3/(D1*D1)); 890 F3 = E1*E1 + D3*D3/(D1*D1) - C; 891 cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n"; 892 if (fabs(F1) < MYEPSILON) { 893 cout << Verbose(3) << "F1 == 0!\n"; 894 cout << Verbose(3) << "Gleichungssystem linear\n"; 895 x[2] = F3/(2.*F2); 896 } else { 897 p = F2/F1; 898 q = p*p - F3/F1; 899 cout << Verbose(3) << "p " << p << "\tq " << q << endl; 900 if (q < 0) { 901 cout << Verbose(3) << "q < 0" << endl; 902 return false; 903 } 904 x[2] = p + sqrt(q); 905 } 906 x[1] = (-D2 * x[2] - D3)/D1; 907 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2]; 908 } 909 switch (flag) { // back-flipping 910 default: 911 case 0: 912 break; 913 case 2: 914 flip(&x1->x[0],&x1->x[1]); 915 flip(&x2->x[0],&x2->x[1]); 916 flip(&y->x[0],&y->x[1]); 917 flip(&x[0],&x[1]); 918 flip(&x1->x[1],&x1->x[2]); 919 flip(&x2->x[1],&x2->x[2]); 920 flip(&y->x[1],&y->x[2]); 921 flip(&x[1],&x[2]); 922 case 1: 923 flip(&x1->x[0],&x1->x[1]); 924 flip(&x2->x[0],&x2->x[1]); 925 flip(&y->x[0],&y->x[1]); 926 //flip(&x[0],&x[1]); 927 flip(&x1->x[1],&x1->x[2]); 928 flip(&x2->x[1],&x2->x[2]); 929 flip(&y->x[1],&y->x[2]); 930 flip(&x[1],&x[2]); 931 break; 932 } 933 // one z component is only determined by its radius (without sign) 934 // thus check eight possible sign flips and determine by checking angle with second vector 935 for (i=0;i<8;i++) { 936 // set sign vector accordingly 937 for (j=2;j>=0;j--) { 938 k = (i & pot(2,j)) << j; 939 cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl; 940 sign[j] = (k == 0) ? 1. : -1.; 941 } 942 cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n"; 943 // apply sign matrix 944 for (j=NDIM;j--;) 945 x[j] *= sign[j]; 946 // calculate angle and check 947 ang = x2->Angle (this); 948 cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t"; 949 if (fabs(ang - cos(beta)) < MYEPSILON) { 950 break; 951 } 952 // unapply sign matrix (is its own inverse) 953 for (j=NDIM;j--;) 954 x[j] *= sign[j]; 955 } 956 return true; 957 }; -
src/vector.hpp
r51c910 r36ec71 8 8 */ 9 9 class Vector { 10 11 10 public: 11 double x[NDIM]; 12 12 13 14 15 13 Vector(); 14 Vector(double x1, double x2, double x3); 15 ~Vector(); 16 16 17 18 19 20 21 22 23 24 17 double Distance(const Vector *y) const; 18 double DistanceSquared(const Vector *y) const; 19 double PeriodicDistance(const Vector *y, const double *cell_size) const; 20 double PeriodicDistanceSquared(const Vector *y, const double *cell_size) const; 21 double ScalarProduct(const Vector *y) const; 22 double Projection(const Vector *y) const; 23 double Norm() const ; 24 double Angle(const Vector *y) const; 25 25 26 void AddVector(const Vector *y); 27 void SubtractVector(const Vector *y); 28 void CopyVector(const Vector *y); 29 void RotateVector(const Vector *y, const double alpha); 30 void VectorProduct(const Vector *y); 31 void ProjectOntoPlane(const Vector *y); 32 void Zero(); 33 void One(double one); 34 void Init(double x1, double x2, double x3); 35 void Normalize(); 36 void Translate(const Vector *x); 37 void Mirror(const Vector *x); 38 void Scale(double **factor); 39 void Scale(double *factor); 40 void Scale(double factor); 41 void MatrixMultiplication(double *M); 42 void InverseMatrixMultiplication(double *M); 43 void KeepPeriodic(ofstream *out, double *matrix); 44 void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors); 45 double CutsPlaneAt(Vector *A, Vector *B, Vector *C); 46 bool GetOneNormalVector(const Vector *x1); 47 bool MakeNormalVector(const Vector *y1); 48 bool MakeNormalVector(const Vector *y1, const Vector *y2); 49 bool MakeNormalVector(const Vector *x1, const Vector *x2, const Vector *x3); 50 bool SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c); 51 bool LSQdistance(Vector **vectors, int dim); 52 void AskPosition(double *cell_size, bool check); 53 bool Output(ofstream *out) const; 26 void AddVector(const Vector *y); 27 void SubtractVector(const Vector *y); 28 void CopyVector(const Vector *y); 29 void RotateVector(const Vector *y, const double alpha); 30 void VectorProduct(const Vector *y); 31 void ProjectOntoPlane(const Vector *y); 32 void Zero(); 33 void One(double one); 34 void Init(double x1, double x2, double x3); 35 void Normalize(); 36 void Translate(const Vector *x); 37 void Mirror(const Vector *x); 38 void Scale(double **factor); 39 void Scale(double *factor); 40 void Scale(double factor); 41 void MatrixMultiplication(double *M); 42 double * InverseMatrix(double *A); 43 void InverseMatrixMultiplication(double *M); 44 void KeepPeriodic(ofstream *out, double *matrix); 45 void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors); 46 double CutsPlaneAt(Vector *A, Vector *B, Vector *C); 47 bool GetOneNormalVector(const Vector *x1); 48 bool MakeNormalVector(const Vector *y1); 49 bool MakeNormalVector(const Vector *y1, const Vector *y2); 50 bool MakeNormalVector(const Vector *x1, const Vector *x2, const Vector *x3); 51 bool SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c); 52 bool LSQdistance(Vector **vectors, int dim); 53 void AskPosition(double *cell_size, bool check); 54 bool Output(ofstream *out) const; 54 55 }; 55 56 -
src/verbose.cpp
r51c910 r36ec71 7 7 ostream& Verbose::print (ostream &ost) const 8 8 { 9 10 11 12 9 for (int i=Verbosity;i--;) 10 ost.put('\t'); 11 //cout << "Verbose(.) called." << endl; 12 return ost; 13 13 }; 14 14
Note:
See TracChangeset
for help on using the changeset viewer.