Changeset e9fa06f

Timestamp:

Dec 18, 2008, 5:21:34 PM (16 years ago)

Author:

Christian Neuen <neuen@…>

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:

7c6712

Parents:

10af0d

Message:

First change to Ball angle as criteria, however, change sucks.

File:

• src/boundary.cpp (modified) (14 diffs)

src/boundary.cpp

@@ -1437 +1437 @@
  * @param mol molecule structure with atoms and bonds
  */
-void
-Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
     atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
@@ -1448 +1447 @@
   Vector dif_a; //Vector from a to candidate
   Vector dif_b; //Vector from b to candidate
-  Vector AngleCheck; // Projection of a difference vector on plane orthogonal on triangle side.
+  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
   Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
 
   double CurrentEpsilon = 0.1;
   double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+  double BallAngle;
   atom *Walker; // variable atom point
 
@@ -1460 +1460 @@
   dif_b.CopyVector(&(b->x));
   dif_b.SubtractVector(&(Candidate->x));
-  AngleCheck.CopyVector(&dif_a);
-  AngleCheck.Scale(-1);
-  AngleCheck.ProjectOntoPlane(Chord);
+  DirectionCheckPoint.CopyVector(&dif_a);
+  DirectionCheckPoint.Scale(-1);
+  DirectionCheckPoint.ProjectOntoPlane(Chord);
 
   SideA = dif_b.Norm();
@@ -1484 +1484 @@
   if (a != Candidate and b != Candidate)
     {
-//      alpha = dif_a.Angle(&dif_b) / 2.;
-//      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-//      SideB = dif_a.Norm();
-//      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-//          alpha); // note this is squared of center line length
-// Those are remains from Freddie. Needed?
-
-//      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-//      Old Chordtest. Replaced by test for Radius of Circumscribing circle
+      //      alpha = dif_a.Angle(&dif_b) / 2.;
+      //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+      //      SideB = dif_a.Norm();
+      //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+      //          alpha); // note this is squared of center line length
+      //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+      // Those are remains from Freddie. Needed?
+
+
 
       Umkreisradius = SideA / 2.0 / sin(alpha);
@@ -1499 +1499 @@
       //cout << SideC / 2.0 / sin(gamma) << endl;
 
-      if (Umkreisradius < RADIUS)
-        {
-
-          sign = AngleCheck.ScalarProduct(d1);
-          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-          // intermediate calculations:
+      if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) //Checking whether ball will at least rest o points.
+        {
+
+          // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
 
           Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
           Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
 
           TempNormal.CopyVector(&dif_a);
@@ -1519 +1515 @@
           TempNormal.Normalize();
           Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-
           TempNormal.Scale(Restradius);
 
@@ -1525 +1520 @@
           Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
 
+          AngleCheck.CopyVector(Chord);
+          AngleCheck.Scale(-0.5);
+          AngleCheck.SubtractVector(&(b->x));
+          AngleCheckReference.CopyVector(&AngleCheck);
+          AngleCheckReference.AddVector(Opt_Mittelpunkt);
+          AngleCheck.AddVector(&Mittelpunkt);
+
+          BallAngle = AngleCheck.Angle(&AngleCheckReference);
+
+          d1->ProjectOntoPlane(&AngleCheckReference);
+          sign = AngleCheck.ScalarProduct(d1);
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
 
 
@@ -1533 +1540 @@
               Opt_Candidate = Candidate;
               Storage[0] = sign;
-              Storage[1] = AngleCheck.Angle(OldNormal);
+              Storage[1] = BallAngle;
               Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
               cout << "Angle is " << Storage[1] << ", Halbraum ist "
@@ -1542 +1549 @@
           else
             {
+              /*
+               * removed due to change in criterium, now checking angle of ball to old normal.
               //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
               //within the ball.
@@ -1550 +1559 @@
 
               if (Distance >RADIUS) // We have no interference and may now check whether the new point is better.
+                */
                 {
                   //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
 
-                  if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon)) || //This will give absolute preference to those in "right-hand" quadrants
-                      (Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
+                  if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon)))  //This will give absolute preference to those in "right-hand" quadrants
+                      //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
                     {
                       cout << "Next better candidate is " << *Candidate << " with ";
                       Opt_Candidate = Candidate;
                       Storage[0] = sign;
-                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Storage[1] = BallAngle;
                       Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
                           << Storage[0] << endl;
-/*
-                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
-                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
-                                    cout << "Restradius ist " << Restradius << endl;
-                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
-                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
-                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
-                                    */
+
 
                     }
@@ -1585 +1579 @@
                     {
                       if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0
-                          && Storage[1] > AngleCheck.Angle(OldNormal)) ||
+                          && Storage[1] > BallAngle) ||
                           (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0
-                          && Storage[1] < AngleCheck.Angle(OldNormal)))
+                          && Storage[1] < BallAngle))
                       //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
                         {
@@ -1593 +1587 @@
                           Opt_Candidate = Candidate;
                           Storage[0] = sign;
-                          Storage[1] = AngleCheck.Angle(OldNormal);
+                          Storage[1] = BallAngle;
                           Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
                           cout << "Angle is " << Storage[1] << ", Halbraum ist "
                               << Storage[0] << endl;
                         }
+
                     }
                 }
-              else
+              /*
+               * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
+               *
+                else
                 {
-                  if (sign>0 && AngleCheck.Angle(OldNormal)>0 && Storage[0]<0)
+                  if (sign>0 && BallAngle>0 && Storage[0]<0)
                     {
                       cout << "Next better candidate is " << *Candidate << " with ";
                       Opt_Candidate = Candidate;
                       Storage[0] = sign;
-                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Storage[1] = BallAngle;
                       Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
                       << Storage[0] << endl;
 
-/*
+//Debugging purposes only
                       cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
-                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
-                                    cout << "Restradius ist " << Restradius << endl;
-                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
-                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
-                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
-                                    */
+                      cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                      cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                      cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                      cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                      cout << "Umkreisradius ist " << Umkreisradius << endl;
+                      cout << "Restradius ist " << Restradius << endl;
+                      cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                      cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                      cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                      cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                      cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+
 
 
@@ -1638 +1636 @@
                     }
                 }
+                */
             }
         }
@@ -1728 +1727 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
+  {
+    Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
+    double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+
+    SideA = dif_b.Norm();
+      SideB = dif_a.Norm();
+      SideC = Chord->Norm();
+      //Chord->Scale(-1);
+
+      alpha = Chord->Angle(&dif_a);
+      beta = M_PI - Chord->Angle(&dif_b);
+      gamma = dif_a.Angle(&dif_b);
+
+      Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+      Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+  }
   cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;