Changeset 110ceb


Ignore:
Timestamp:
Jun 12, 2008, 7:14:46 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
6c5812
Parents:
bc84e47
Message:

VolumeOfConvexEnvelope: Works!

VolumeOfConvexEnvelope has been analysed into various smaller functions and approach is working.
two new files: boundary.?pp
various new functions:
class Tesselation with AddPoint(), TesselateOnBoundary() and GuessStartingTriangle() does the actual tesselation
CreateClustersinWater() will create the repetition of the cluster with correct spacing (unfinished).
GetDiametersOfCluster() calculate the greatest diameter in projection per axis
GetBoundaryPoints() gets the boundary on the convex envelope by projection for a molecular cluster
GetCommonEndpoint() finds the endpoint two lines are sharing

Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rbc84e47 r110ceb  
    1 SOURCE = atom.cpp bond.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
    2 HEADER = defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
     1SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
     2HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp stackclass.hpp vector.hpp
    33
    44bin_PROGRAMS = molecuilder joiner analyzer

  • src/builder.cpp

    rbc84e47 r110ceb  
    5252#include "helpers.hpp"
    5353#include "molecules.hpp"
     54#include "boundary.hpp"
    5455
    5556/********************************************** Submenu routine **************************************/
    … …  
    566567    case 'e':
    567568        cout << Verbose(0) << "Evaluating volume of the convex envelope.";
    568         mol->VolumeOfConvexEnvelope((ofstream *)&cout, configuration->GetIsAngstroem());
     569        VolumeOfConvexEnvelope((ofstream *)&cout, configuration, mol);
    569570        break;
    570571  }
    … …  
    733734  string line;
    734735  atom *first;
     736  double tmp;
    735737  int ExitFlag = 0;
    736738  int j;
    … …  
    762764            cout << "\t-m\tAlign in PAS with greatest EV along z axis." << endl;
    763765            cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
     766            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    764767            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
    765768            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    … …  
    967970            case 'm':
    968971              ExitFlag = 1;
    969               cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
    970               mol->PrincipalAxisSystem((ofstream *)&cout, true);
     972              j = atoi(argv[argptr++]);
     973              if ((j<0) || (j>1)) {
     974                cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
     975                j = 0;
     976              }
     977              if (j)
     978                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
     979              else
     980                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
     981              mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
     982              break;
     983            case 'o':
     984              ExitFlag = 1;
     985              cout << Verbose(0) << "Evaluating volume of the convex envelope.";
     986//              tmp = atof(argv[argptr++]);
     987//              if (tmp < 1.0) {
     988//                cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
     989//                tmp = 1.3;
     990//              }
     991              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, mol);
    971992              break;
    972993            default:   // no match? Step on
    973               argptr++;
     994              if (argv[argptr][0] != '-') // if it started with a '-' we've already made a step!
     995                argptr++;
    974996              break;
    975997          }
    … …  
    12451267 
    12461268  // save element data base
    1247   if (periode->StorePeriodentafel()) //ElementsFileName
     1269  if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
    12481270    cout << Verbose(0) << "Saving of elements.db successful." << endl;
    12491271  else

  • src/defs.hpp

    rbc84e47 r110ceb  
    4343
    4444// various standard filenames
    45 #define DEFAULTCONFIG "main_pcp_linux"
    46 #define KEYSETFILE "KeySets.dat"
    47 #define ADJACENCYFILE "Adjacency.dat"
    48 #define TEFACTORSFILE "TE-Factors.dat"
    49 #define FORCESFILE "Forces-Factors.dat"
    50 #define ORDERATSITEFILE "OrderAtSite.dat"
    51 #define ENERGYPERFRAGMENT "EnergyPerFragment"
    52 #define FRAGMENTPREFIX "BondFragment"
    53 #define STANDARDCONFIG "unknown.conf"
    54 #define STANDARDELEMENTSDB "elements.db"
    55 #define STANDARDVALENCEDB "valence.db"
    56 #define STANDARDORBITALDB "orbitals.db"
    57 #define STANDARDHBONDDISTANCEDB "Hbonddistance.db"
    58 #define STANDARDHBONDANGLEDB "Hbondangle.db"
     45#define DEFAULTCONFIG "main_pcp_linux"    //!< default filename of config file
     46#define CONVEXENVELOPE "ConvexEnvelope.dat"    //!< default filename of convex envelope tecplot data file
     47#define KEYSETFILE "KeySets.dat"    //!< default filename of BOSSANOVA key sets file
     48#define ADJACENCYFILE "Adjacency.dat"    //!< default filename of BOSSANOVA adjacancy file
     49#define TEFACTORSFILE "TE-Factors.dat"    //!< default filename of BOSSANOVA total energy factors file
     50#define FORCESFILE "Forces-Factors.dat"    //!< default filename of BOSSANOVA force factors file
     51#define ORDERATSITEFILE "OrderAtSite.dat"    //!< default filename of BOSSANOVA Bond Order at each atom file
     52#define ENERGYPERFRAGMENT "EnergyPerFragment"    //!< default filename of BOSSANOVA Energy contribution Per Fragment file
     53#define FRAGMENTPREFIX "BondFragment"    //!< default filename prefix of BOSSANOVA fragment config and directories
     54#define STANDARDCONFIG "unknown.conf"    //!< default filename of standard config file
     55#define STANDARDELEMENTSDB "elements.db"    //!< default filename of elements data base with masses, Z, VanDerWaals radii, ...
     56#define STANDARDVALENCEDB "valence.db"    //!< default filename of valence number per element database
     57#define STANDARDORBITALDB "orbitals.db"    //!< default filename of orbitals per element database
     58#define STANDARDHBONDDISTANCEDB "Hbonddistance.db"    //!< default filename of typial bond distance to hydrogen database
     59#define STANDARDHBONDANGLEDB "Hbondangle.db"    //!< default filename of typial bond angle to hydrogen database
     60
     61// some values
     62#define SOLVENTDENSITY_A 0.6022142
     63#define SOLVENTDENSITY_a0 0.089238936
     64
    5965
    6066#define UPDATECOUNT 10  //!< update ten sites per BOSSANOVA interval

  • src/molecules.cpp

    rbc84e47 r110ceb  
    13701370  return No;
    13711371};
    1372 
    1373 /** Determines the volume of a cluster.
    1374  * Determines first the convex envelope, then tesselates it and calculates its volume.
    1375  * \param *out output stream for debugging
    1376  * \param IsAngstroem for the units on output
    1377  */
    1378 double molecule::VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem)
    1379 {
    1380         atom *Walker = NULL;
    1381        
    1382         // 1. calculate center of gravity
    1383         *out << endl;
    1384         vector *CenterOfGravity = DetermineCenterOfGravity(out);
    1385        
    1386         // 2. translate all points into CoG
    1387         *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
    1388         Walker=start;
    1389         while (Walker->next != end) {
    1390           Walker = Walker->next;
    1391           Walker->x.Translate(CenterOfGravity);
    1392         }
    1393         // 2. calculate distance of each atom and sort into hash table
    1394 //      map<double, int> DistanceSet;
    1395 //     
    1396 //      Walker = start;
    1397 //      while (Walker->next != end) {
    1398 //              Walker = Walker->next;
    1399 //              DistanceSet.insert( pair<double, int>(ptr->x.distance(CenterOfGravity), ptr->nr) );
    1400 //      }
    1401        
    1402         // 3. Find all points on the boundary
    1403   *out << Verbose(1) << "Finding all boundary points." << endl;
    1404         Boundaries BoundaryPoints[NDIM];        // first is alpha, second is (r, nr)
    1405         BoundariesTestPair BoundaryTestPair;
    1406         vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
    1407         double radius, angle;
    1408         // 3a. Go through every axis
    1409         for (int axis=0; axis<NDIM; axis++)  {
    1410                 AxisVector.Zero();
    1411     AngleReferenceVector.Zero();
    1412     AngleReferenceNormalVector.Zero();
    1413                 AxisVector.x[axis] = 1.;
    1414                 AngleReferenceVector.x[(axis+1)%NDIM] = 1.;
    1415                 AngleReferenceNormalVector.x[(axis+2)%NDIM] = 1.;
    1416 //              *out << Verbose(1) << "Axisvector is ";
    1417 //              AxisVector.Output(out);
    1418 //              *out << " and AngleReferenceVector is ";
    1419 //              AngleReferenceVector.Output(out);
    1420 //              *out << "." << endl;
    1421 //    *out << " and AngleReferenceNormalVector is ";
    1422 //    AngleReferenceNormalVector.Output(out);
    1423 //    *out << "." << endl;
    1424                 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    1425                 Walker = start;
    1426                 while (Walker->next != end) {
    1427                         Walker = Walker->next;
    1428                         vector ProjectedVector;
    1429                         ProjectedVector.CopyVector(&Walker->x);
    1430                         ProjectedVector.ProjectOntoPlane(&AxisVector);
    1431                   // correct for negative side
    1432                   //if (Projection(y) < 0)
    1433                     //angle = 2.*M_PI - angle;
    1434                         radius = ProjectedVector.Norm();
    1435                         if (fabs(radius) > MYEPSILON)
    1436                           angle = ProjectedVector.Angle(&AngleReferenceVector);
    1437                         else
    1438                           angle = 0.;  // otherwise it's a vector in Axis Direction and unimportant for boundary issues
    1439                          
    1440                         //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
    1441                         if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
    1442                           angle = 2.*M_PI - angle;
    1443                         }
    1444                         *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
    1445                         ProjectedVector.Output(out);
    1446                         *out << endl;
    1447                         BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (angle, DistanceNrPair (radius, Walker) ) );
    1448                         if (BoundaryTestPair.second) { // successfully inserted
    1449                         } else { // same point exists, check first r, then distance of original vectors to center of gravity
    1450                                 *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
    1451                                 *out << Verbose(2) << "Present vector: ";
    1452                                 BoundaryTestPair.first->second.second->x.Output(out);
    1453                                 *out << endl;
    1454                                 *out << Verbose(2) << "New vector: ";
    1455                                 Walker->x.Output(out);
    1456                                 *out << endl;
    1457                                 double tmp = ProjectedVector.Norm();
    1458                                 if (tmp > BoundaryTestPair.first->second.first) {
    1459                                         BoundaryTestPair.first->second.first = tmp;
    1460                                         BoundaryTestPair.first->second.second = Walker;
    1461                                         *out << Verbose(2) << "Keeping new vector." << endl;
    1462                                 } else if (tmp == BoundaryTestPair.first->second.first) {
    1463                                         if (BoundaryTestPair.first->second.second->x.ScalarProduct(&BoundaryTestPair.first->second.second->x) < Walker->x.ScalarProduct(&Walker->x)) { // Norm() does a sqrt, which makes it a lot slower
    1464                                                 BoundaryTestPair.first->second.second = Walker;
    1465                                                 *out << Verbose(2) << "Keeping new vector." << endl;
    1466                                         } else {
    1467                                                 *out << Verbose(2) << "Keeping present vector." << endl;
    1468                                         }
    1469                                 } else {
    1470                                                 *out << Verbose(2) << "Keeping present vector." << endl;
    1471                                 }
    1472                         }
    1473                 }
    1474                 // printing all inserted for debugging
    1475                 {
    1476       *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    1477       int i=0;
    1478       for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1479         if (runner != BoundaryPoints[axis].begin())
    1480           *out << ", " << i << ": " << *runner->second.second;
    1481         else
    1482           *out << i << ": " << *runner->second.second;
    1483         i++;
    1484       }
    1485       *out << endl;
    1486                 }
    1487     // 3c. throw out points whose distance is less than the mean of left and right neighbours
    1488                 bool flag = false;
    1489                 do { // do as long as we still throw one out per round
    1490                   *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
    1491                         flag = false;
    1492                         Boundaries::iterator left = BoundaryPoints[axis].end();
    1493                         Boundaries::iterator right = BoundaryPoints[axis].end();
    1494                         for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1495                                 // set neighbours correctly
    1496                                 if (runner == BoundaryPoints[axis].begin()) {
    1497                                         left = BoundaryPoints[axis].end();
    1498                                 }       else {
    1499                                         left = runner;
    1500                                 }
    1501         left--;
    1502         right = runner;
    1503         right++;
    1504                                 if (right == BoundaryPoints[axis].end()) {
    1505                                         right = BoundaryPoints[axis].begin();
    1506                                 }
    1507                                 // check distance
    1508                                
    1509                                 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
    1510                                 {
    1511           vector SideA, SideB, SideC, SideH;
    1512           SideA.CopyVector(&left->second.second->x);
    1513           SideA.ProjectOntoPlane(&AxisVector);
    1514 //          *out << "SideA: ";
    1515 //          SideA.Output(out);
    1516 //          *out << endl;
    1517          
    1518           SideB.CopyVector(&right->second.second->x);
    1519           SideB.ProjectOntoPlane(&AxisVector);
    1520 //          *out << "SideB: ";
    1521 //          SideB.Output(out);
    1522 //          *out << endl;
    1523          
    1524           SideC.CopyVector(&left->second.second->x);
    1525           SideC.SubtractVector(&right->second.second->x);
    1526           SideC.ProjectOntoPlane(&AxisVector);
    1527 //          *out << "SideC: ";
    1528 //          SideC.Output(out);
    1529 //          *out << endl;
    1530 
    1531           SideH.CopyVector(&runner->second.second->x);
    1532           SideH.ProjectOntoPlane(&AxisVector);
    1533 //          *out << "SideH: ";
    1534 //          SideH.Output(out);
    1535 //          *out << endl;
    1536          
    1537           // calculate each length
    1538           double a = SideA.Norm();
    1539           //double b = SideB.Norm();
    1540           //double c = SideC.Norm();
    1541           double h = SideH.Norm();
    1542           // calculate the angles
    1543           double alpha = SideA.Angle(&SideH);
    1544                                 double beta = SideA.Angle(&SideC);
    1545                                 double gamma = SideB.Angle(&SideH);
    1546                                 double delta = SideC.Angle(&SideH);
    1547                                 double MinDistance = a * sin(beta)/(sin(delta)) * (((alpha < M_PI/2.) || (gamma < M_PI/2.)) ? 1. : -1.);
    1548 //                              *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
    1549                                 *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
    1550                                 if ((fabs(h/fabs(h) - MinDistance/fabs(MinDistance)) < MYEPSILON) && (h <  MinDistance)) {
    1551                                         // throw out point
    1552                                   *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
    1553                                         BoundaryPoints[axis].erase(runner);
    1554                                         flag = true;
    1555                                 }
    1556                                 }
    1557                         }
    1558                 } while (flag);
    1559         }       
    1560         // 3d. put into boundary set only those points appearing in each of the NDIM sets
    1561         int *AtomList = new int[AtomCount];
    1562         for (int j=0; j<AtomCount; j++) // reset list
    1563                 AtomList[j] = 0;
    1564         for (int axis=0; axis<NDIM; axis++)  {  // increase whether it's present
    1565                 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1566                         AtomList[runner->second.second->nr]++;
    1567                 }
    1568         }
    1569 
    1570         // 3e. Points no more in the total list, have to be thrown out of each axis lists too!
    1571         for (int axis=0; axis<NDIM; axis++)  {
    1572                 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1573                         if (AtomList[runner->second.second->nr] < 1) {
    1574                           *out << Verbose(1) << "Throwing out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
    1575                                 BoundaryPoints[axis].erase(runner);
    1576                         }
    1577                 }
    1578         }
    1579        
    1580         // listing boundary points
    1581         *out << Verbose(1) << "The following atoms are on the boundary:";
    1582   int BoundaryPointCount = 0;
    1583   Walker = start;
    1584   while (Walker->next != end) {
    1585     Walker = Walker->next;
    1586     if (AtomList[Walker->nr] > 0) {
    1587       BoundaryPointCount ++;
    1588       *out << " " << Walker->Name;
    1589     }
    1590   }
    1591   *out << endl;
    1592 
    1593         *out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
    1594         // now we have the whole set of edge points in the BoundaryList
    1595        
    1596         // 4. Create each possible line, number them and put them into a list (3 * AtomCount possible)
    1597         LinesMap AllTriangleLines;
    1598        
    1599         // then create each line
    1600   *out << Verbose(1) << "Creating lines between adjacent boundary points." << endl;
    1601   atom *Runner = NULL;
    1602   int LineCount = 0;
    1603   Walker = start;
    1604   while (Walker->next != end) {
    1605     Walker = Walker->next;
    1606     if (AtomList[Walker->nr] > 0) { // boundary point!
    1607       *out << Verbose(1) << "Current Walker is " << *Walker << "." << endl;
    1608       // make a list of all neighbours
    1609       DistanceMap Distances;
    1610       double distance = 0;
    1611       Runner = start;
    1612       while (Runner->next != end) {
    1613         Runner = Runner->next;
    1614         if ((AtomList[Runner->nr] > 0) && (Runner != Walker)) { // don't mess with ourselves!
    1615           distance = Walker->x.Distance(&Runner->x);    // note this is squared distance
    1616           Distances.insert ( DistanceNrPair( distance, Runner) );
    1617           *out << Verbose(2) << "Inserted distance " << distance << " to " << *Runner << " successfully." << endl;
    1618         }
    1619       } // end of distance filling while loop
    1620       // take 3NNs to draw a line in between
    1621       int Counter = 0;  // counts up to three inserted lines
    1622       LinesTestPair LinesTest;
    1623       for (DistanceMap::iterator spanner = Distances.begin(); spanner != Distances.end(); spanner++) {
    1624         *out << Verbose(1) << "Current distance to insert is " << spanner->first << " to " << spanner->second+1 << "." << endl;
    1625         LinesTest = AllTriangleLines.insert ( LinesPair( (Walker->nr*AtomCount)+spanner->second->nr, LineCount) );
    1626         if (!LinesTest.second) {
    1627           *out << Verbose(2) << "Insertion of line between " << Walker->nr+1 << " and " << spanner->second->nr+1 << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
    1628         } else {
    1629           LineCount++;
    1630           Counter++;
    1631           *out << Verbose(2) << "Insertion of line between " << (LinesTest.first->first/AtomCount)+1 << " and " << (LinesTest.first->first%AtomCount)+1 << " successful, inserting mirrored line and increased counter to " << Counter << "." << endl;
    1632           LinesTest = AllTriangleLines.insert ( LinesPair( Walker->nr+(spanner->second->nr*AtomCount), LineCount) );
    1633           if (!LinesTest.second) {
    1634             *out << Verbose(2) << "Insertion of mirror line between " << Walker->nr << " and " << spanner->second->nr << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
    1635             *out << Verbose(0) << "SERIOUS ERROR!!!" << endl;
    1636           }
    1637         }
    1638         if (Counter >= 3) { // stop after three lines (if not run out of possible NNs already)
    1639           *out << Verbose(2) << "Ok, three lines inserted from current walker " << *Walker << ", stopping." << endl;
    1640           break;
    1641         }
    1642       }
    1643     } // end of if boundary pointy
    1644   } // end of loop through all atoms
    1645 
    1646         // listing created lines
    1647         *out << Verbose(2) << "The following lines were created:";
    1648         for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++)
    1649           if ((liner->first/AtomCount) < (liner->first%AtomCount))
    1650             *out << " " << (liner->first/AtomCount)+1 << "<->" << (liner->first%AtomCount)+1;
    1651         *out << endl;
    1652         *out << Verbose(2) << "I created " << LineCount << " edges." << endl;
    1653        
    1654         // 5. Create every possible triangle from the lines (numbering of each line must be i < j < k)
    1655         struct triangles {
    1656                 atom *x[3];
    1657                 int nr;
    1658         } TriangleList[2 * 3 * AtomCount];      // each line can be part in at most 2 triangles
    1659 
    1660   *out << Verbose(1) << "Creating triangles out of these lines." << endl;
    1661         int TriangleCount = 0;
    1662   LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
    1663   int endpoints[3];
    1664   for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++) { // go through every line
    1665                 LinetoAtomMap LinetoAtom;
    1666                 // we have two points, check the lines at either end, whether they share a common endpoint
    1667     *out << Verbose(1) << "Examining line between " << (liner->first/AtomCount)+1 << " and " << (liner->first%AtomCount)+1 << "." << endl;
    1668     int enden[3][2];
    1669     enden[0][0] = (liner->first/AtomCount);
    1670     enden[0][1] = (liner->first%AtomCount);
    1671     if (enden[0][0] < enden[0][1]) {
    1672                 for (int endpoint=0;endpoint<2;endpoint++) {
    1673                         endpoints[0] = enden[0][endpoint];
    1674                         *out << Verbose(2) << "Current first endpoint is " << endpoints[0]+1 << "." << endl;
    1675                        
    1676                         LinesMap::iterator LineRangeStart = AllTriangleLines.lower_bound(endpoints[0]*AtomCount);
    1677                         for (LinesMap::iterator coach = LineRangeStart; (coach->first/AtomCount) == endpoints[0]; coach++) { // look at all line's other endpoint
    1678                     enden[1][0] = (coach->first/AtomCount);
    1679                     enden[1][1] = (coach->first%AtomCount);
    1680           //*out << Verbose(3) << "Line is " << coach->first << " with nr. " << coach->second << ": " << enden[1][0+1] << "<->" << enden[1][1]+1 << "." << endl;
    1681                           endpoints[1] = (enden[1][0] != endpoints[0] ) ? enden[1][0] : enden[1][1];
    1682                           if ((endpoints[1] > endpoints[0]) && (endpoints[1] != enden[0][(endpoint+1)%2])) {  // don't go back the very same line!
    1683                           *out << Verbose(3) << "Current second endpoint is " << endpoints[1]+1 << "." << endl; 
    1684                          
    1685               LinesMap::iterator LineRangeStart2 = AllTriangleLines.lower_bound(endpoints[1]*AtomCount);
    1686               for (LinesMap::iterator coacher = LineRangeStart2; (coacher->first/AtomCount) == endpoints[1]; coacher++) { // look at all line's other endpoint
    1687                 enden[2][0] = (coacher->first/AtomCount);
    1688                 enden[2][1] = (coacher->first%AtomCount);
    1689               //*out << Verbose(3) << "Line is " << coacher->first << " with nr. " << coacher->second << ": " << enden[2][0]+1 << "<->" << enden[2][1]+1 << "." << endl;
    1690                 endpoints[2] = (enden[2][0] != endpoints[1]) ? enden[2][0] : enden[2][1];
    1691                 if ((endpoints[2] > endpoints[1]) && (endpoints[2] != endpoints[0])) {  // don't go back the very same line!
    1692                 *out << Verbose(4) << "Current third endpoint is " << endpoints[2]+1 << "." << endl;
    1693                
    1694                 if (endpoints[2] == enden[0][(endpoint+1)%2]) { // jo, closing triangle
    1695                   *out << Verbose(3) << "MATCH: Triangle is ";
    1696                   for (int k=0;k<3;k++) {
    1697                     *out << endpoints[k]+1 << "\t";
    1698                     TriangleList[TriangleCount].x[k] = FindAtom(endpoints[k]);
    1699                   }
    1700                   *out << endl;
    1701                   TriangleList[TriangleCount].nr = TriangleCount;
    1702                   TriangleCount++;
    1703                 } else {
    1704                   *out << Verbose(3) << "No match!" << endl;
    1705                 }
    1706               }
    1707               }
    1708                         }
    1709                         }
    1710                 }
    1711     }
    1712         }
    1713   // listing created lines
    1714   *out << Verbose(2) << "The following triangles were created:";
    1715   for (int i=0;i<TriangleCount;i++) {
    1716     *out << " " << TriangleList[i].x[0]->Name << "<->" << TriangleList[i].x[1]->Name << "<->" << TriangleList[i].x[2]->Name;
    1717   }
    1718   *out << endl;
    1719         *out << Verbose(2) << "I created " << TriangleCount << " triangles." << endl;
    1720 
    1721         // 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
    1722   *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
    1723         double volume = 0.;
    1724         double PyramidVolume = 0.;
    1725         double G,h;
    1726         vector x,y;
    1727         double a,b,c;
    1728         for (int i=0; i<TriangleCount; i++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
    1729                 x.CopyVector(&TriangleList[i].x[0]->x);
    1730                 x.SubtractVector(&TriangleList[i].x[1]->x);
    1731                 y.CopyVector(&TriangleList[i].x[0]->x);
    1732                 y.SubtractVector(&TriangleList[i].x[2]->x);
    1733                 a = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[1]->x));
    1734                 b = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[2]->x));
    1735                 c = sqrt(TriangleList[i].x[2]->x.Distance(&TriangleList[i].x[1]->x));
    1736                 G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
    1737                 x.MakeNormalVector(&TriangleList[i].x[0]->x, &TriangleList[i].x[1]->x, &TriangleList[i].x[2]->x);
    1738                 x.Scale(TriangleList[i].x[1]->x.Projection(&x));
    1739                 h = x.Norm(); // distance of CoG to triangle
    1740                 PyramidVolume = (1./3.) * G * h;                // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
    1741                 *out << Verbose(2) << "Area of triangle is " << G << "A^2, height is " << h << " and the volume is " << PyramidVolume << "." << endl;
    1742                 volume += PyramidVolume;
    1743         }
    1744         *out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
    1745 
    1746   // 7. translate all points back from CoG
    1747   *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
    1748         CenterOfGravity->Scale(-1);
    1749   Walker=start;
    1750   while (Walker->next != end) {
    1751     Walker = Walker->next;
    1752     Walker->x.Translate(CenterOfGravity);
    1753   }
    1754 
    1755   // free reference lists
    1756         delete[](AtomList);
    1757 
    1758         return volume;
    1759 };
    1760 
    17611372/** Returns Shading as a char string.
    17621373 * \param color the Shading

  • src/molecules.hpp

    rbc84e47 r110ceb  
    4444#define KeySetTestPair pair<KeySet::iterator, bool>
    4545#define GraphTestPair pair<Graph::iterator, bool>
    46 
    47 #define DistanceNrPair pair < double, atom* >
    48 #define DistanceMap multimap < double, atom* >
    49 #define DistanceTestPair pair < DistanceMap::iterator, bool>
    50 
    51 #define Boundaries map <double, DistanceNrPair >
    52 #define BoundariesPair pair<double, DistanceNrPair >
    53 #define BoundariesTestPair pair< Boundaries::iterator, bool>
    54 
    55 #define LinesMap map <int, int>
    56 #define LinesPair pair <int, int>
    57 #define LinesTestPair pair < LinesMap::iterator, bool>
    58 
    59 #define LinetoAtomMap map < int, atom * >
    60 #define LinetoAtomPair pair < int, atom * >
    61 #define LinetoAtomTestPair pair < LinetoAtomMap::iterator, bool>
    6246
    6347struct KeyCompare

  • src/vector.cpp

    rbc84e47 r110ceb  
    418418    return false;
    419419  }
    420   cout << Verbose(4) << "relative, first plane coordinates:";
    421   x1.Output((ofstream *)&cout);
    422   cout << endl;
    423   cout << Verbose(4) << "second plane coordinates:";
    424   x2.Output((ofstream *)&cout);
    425   cout << endl;
     420//  cout << Verbose(4) << "relative, first plane coordinates:";
     421//  x1.Output((ofstream *)&cout);
     422//  cout << endl;
     423//  cout << Verbose(4) << "second plane coordinates:";
     424//  x2.Output((ofstream *)&cout);
     425//  cout << endl;
    426426
    427427  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    … …  
    452452    return false;
    453453  }
    454   cout << Verbose(4) << "relative, first plane coordinates:";
    455   x1.Output((ofstream *)&cout);
    456   cout << endl;
    457   cout << Verbose(4) << "second plane coordinates:";
    458   x2.Output((ofstream *)&cout);
    459   cout << endl;
     454//  cout << Verbose(4) << "relative, first plane coordinates:";
     455//  x1.Output((ofstream *)&cout);
     456//  cout << endl;
     457//  cout << Verbose(4) << "second plane coordinates:";
     458//  x2.Output((ofstream *)&cout);
     459//  cout << endl;
    460460
    461461  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    … …  
    529529      return false;
    530530  }
     531};
     532
     533/** Determines paramter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
     534 * \param *A first plane vector
     535 * \param *B second plane vector
     536 * \param *C third plane vector
     537 * \return scaling parameter for this vector
     538 */
     539double vector::CutsPlaneAt(vector *A, vector *B, vector *C)
     540{
     541//  cout << Verbose(3) << "For comparison: ";
     542//  cout << "A " << A->Projection(this) << "\t";
     543//  cout << "B " << B->Projection(this) << "\t";
     544//  cout << "C " << C->Projection(this) << "\t";
     545//  cout << endl;
     546  return A->Projection(this);
    531547};
    532548

  • src/vector.hpp

    rbc84e47 r110ceb  
    3939  void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
    4040 
     41  double CutsPlaneAt(vector *A, vector *B, vector *C);
    4142  bool GetOneNormalVector(const vector *x1);
    4243  bool MakeNormalVector(const vector *y1);
Note: See TracChangeset for help on using the changeset viewer.