Changeset 6c438f for src/Box.cpp

Timestamp:

Aug 28, 2010, 3:21:11 AM (15 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:

a6e6b5, f8982c

Parents:

2ad482 (diff), fd4905 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Frederik Heber <heber@…> (08/28/10 03:17:48)

git-committer:

Frederik Heber <heber@…> (08/28/10 03:21:11)

Message:

Merge branch 'StructureRefactoring' into Shapes

Conflicts:

src/Box.cpp
src/Box.hpp
src/Descriptors/AtomShapeDescriptor.cpp
src/Descriptors/AtomShapeDescriptor.hpp
src/Descriptors/AtomShapeDescriptor_impl.hpp
src/LinearAlgebra/Line.cpp
src/LinearAlgebra/Line.hpp
src/LinearAlgebra/Matrix.cpp
src/LinearAlgebra/Matrix.hpp
src/Makefile.am
src/Shapes/BaseShapes.cpp
src/Shapes/BaseShapes_impl.hpp
src/Shapes/Shape.cpp
src/Shapes/Shape.hpp
src/Shapes/ShapeOps_impl.hpp
src/Shapes/Shape_impl.hpp
src/unittests/ShapeUnittest.cpp

• AtomShapeDescriptor had been implemented doubly.
• LineSegment* moved to LinearAlgebra/
• BaseShapes:: Cuboid is from (0,0,0) and not from (-1,-1,-1) to (1,1,1)
  • had to fix ShapeUnittest
• some overlap in Shapes due to getHomogeneousPointsOnSurface added in stable.
• BUGFIX: Matrix::transpose() did not allocate MatrixContent::content.
• merged BaseShapesUnitTest into ShapeUnitTest.

File:

• src/Box.cpp (modified) (9 diffs)

src/Box.cpp

@@ +1 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
 /*
  * Box.cpp
@@ -6 +13 @@
  */
 
-//#include "Helpers/MemDebug.hpp"
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "Helpers/MemDebug.hpp"
 
 #include "Box.hpp"
 
 #include <cmath>
-
-#include "Matrix.hpp"
-#include "vector.hpp"
+#include <iostream>
+#include <cstdlib>
+
+#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Plane.hpp"
 #include "Shapes/BaseShapes.hpp"
 #include "Shapes/ShapeOps.hpp"
 
 #include "Helpers/Assert.hpp"
+
+using namespace std;
 
 Box::Box()
@@ -25 +42 @@
   Minv = new Matrix();
   Minv->one();
+  conditions.resize(3);
+  conditions[0] = conditions[1] = conditions[2] = Wrap;
 }
 
@@ -30 +49 @@
   M=new Matrix(*src.M);
   Minv = new Matrix(*src.Minv);
+  conditions = src.conditions;
 }
 
@@ -62 +82 @@
   Vector helper = translateOut(point);
   for(int i=NDIM;i--;){
-    double intpart,fracpart;
-    fracpart = modf(helper.at(i),&intpart);
-    if(fracpart<0.)
-      fracpart+=1.;
-    helper.at(i)=fracpart;
+
+    switch(conditions[i]){
+    case Wrap:
+      helper.at(i)=fmod(helper.at(i),1);
+      helper.at(i)+=(helper.at(i)>=0)?0:1;
+      break;
+    case Bounce:
+      {
+        // there probably is a better way to handle this...
+        // all the fabs and fmod modf probably makes it very slow
+        double intpart,fracpart;
+        fracpart = modf(fabs(helper.at(i)),&intpart);
+        helper.at(i) = fabs(fracpart-fmod(intpart,2));
+      }
+      break;
+    case Ignore:
+      break;
+    default:
+      ASSERT(0,"No default case for this");
+    }
+
   }
   return translateIn(helper);
@@ -74 +110 @@
 {
   bool result = true;
-  Vector tester = (*Minv) * point;
+  Vector tester = translateOut(point);
 
   for(int i=0;i<NDIM;i++)
-    result = result && ((tester[i] >= -MYEPSILON) && ((tester[i] - 1.) < MYEPSILON));
+    result = result &&
+              ((conditions[i] == Ignore) ||
+               ((tester[i] >= -MYEPSILON) &&
+               ((tester[i] - 1.) < MYEPSILON)));
 
   return result;
@@ -84 +123 @@
 
 VECTORSET(std::list) Box::explode(const Vector &point,int n) const{
+  ASSERT(isInside(point),"Exploded point not inside Box");
   VECTORSET(std::list) res;
 
-  // translate the Vector into each of the 27 neighbourhoods
-
-  // first create all translation Vectors
-  // there are (n*2+1)^3 such vectors
-  int max_dim = (n*2+1);
-  int max_dim2 = max_dim*max_dim;
-  int max = max_dim2*max_dim;
-  // only one loop to avoid unneccessary jumps
-  for(int i = 0;i<max;++i){
-    // get all coordinates for this iteration
-    int n1 = (i%max_dim)-n;
-    int n2 = ((i/max_dim)%max_dim)-n;
-    int n3 = ((i/max_dim2))-n;
-    Vector translation = translateIn(Vector(n1,n2,n3));
-    res.push_back(translation);
-  }
-  // translate all the translation vector by the offset defined by point
-  res.translate(point);
+  Vector translater = translateOut(point);
+  Vector mask; // contains the ignored coordinates
+
+  // count the number of coordinates we need to do
+  int dims = 0; // number of dimensions that are not ignored
+  vector<int> coords;
+  vector<int> index;
+  for(int i=0;i<NDIM;++i){
+    if(conditions[i]==Ignore){
+      mask[i]=translater[i];
+      continue;
+    }
+    coords.push_back(i);
+    index.push_back(-n);
+    dims++;
+  } // there are max vectors in total we need to create
+
+  if(!dims){
+    // all boundaries are ignored
+    res.push_back(point);
+    return res;
+  }
+
+  bool done = false;
+  while(!done){
+    // create this vector
+    Vector helper;
+    for(int i=0;i<dims;++i){
+      switch(conditions[coords[i]]){
+        case Wrap:
+          helper[coords[i]] = index[i]+translater[coords[i]];
+          break;
+        case Bounce:
+          {
+            // Bouncing the coordinate x produces the series:
+            // 0 -> x
+            // 1 -> 2-x
+            // 2 -> 2+x
+            // 3 -> 4-x
+            // 4 -> 4+x
+            // the first number is the next bigger even number (n+n%2)
+            // the next number is the value with alternating sign (x-2*(n%2)*x)
+            // the negative numbers produce the same sequence reversed and shifted
+            int n = abs(index[i]) + ((index[i]<0)?-1:0);
+            int sign = (index[i]<0)?-1:+1;
+            int even = n%2;
+            helper[coords[i]]=n+even+translater[coords[i]]-2*even*translater[coords[i]];
+            helper[coords[i]]*=sign;
+          }
+          break;
+        case Ignore:
+          ASSERT(0,"Ignored coordinate handled in generation loop");
+        default:
+          ASSERT(0,"No default case for this switch-case");
+      }
+
+    }
+    // add back all ignored coordinates (not handled in above loop)
+    helper+=mask;
+    res.push_back(translateIn(helper));
+    // set the new indexes
+    int pos=0;
+    ++index[pos];
+    while(index[pos]>n){
+      index[pos++]=-n;
+      if(pos>=dims) { // it's trying to increase one beyond array... all vectors generated
+        done = true;
+        break;
+      }
+      ++index[pos];
+    }
+  }
   return res;
 }
 
 VECTORSET(std::list) Box::explode(const Vector &point) const{
-  VECTORSET(std::list) res;
-
-  // translate the Vector into each of the 27 neighbourhoods
-
-  // first create all 27 translation Vectors
-  // these loops depend on fixed parameters and can easily be expanded
-  // by the compiler to allow code without jumps
-  for(int n1 = -1;n1<=1;++n1){
-    for(int n2 = -1;n2<=1;++n2){
-      for(int n3 = -1;n3<=1;++n3){
-        // get all coordinates for this iteration
-        Vector translation = translateIn(Vector(n1,n2,n3));
-        res.push_back(translation);
-      }
-    }
-  }
-  // translate all the translation vector by the offset defined by point
-  res.translate(point);
-  return res;
+  ASSERT(isInside(point),"Exploded point not inside Box");
+  return explode(point,1);
 }
 
 double Box::periodicDistanceSquared(const Vector &point1,const Vector &point2) const{
-  VECTORSET(std::list) expansion = explode(point1);
-  double res = expansion.minDistSquared(point2);
+  Vector helper1 = WrapPeriodically(point1);
+  Vector helper2 = WrapPeriodically(point2);
+  VECTORSET(std::list) expansion = explode(helper1);
+  double res = expansion.minDistSquared(helper2);
   return res;
 }
@@ -143 +222 @@
 Shape Box::getShape() const{
   return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
-
+}
+
+const Box::Conditions_t Box::getConditions(){
+  return conditions;
+}
+
+void Box::setCondition(int i,Box::BoundaryCondition_t condition){
+  conditions[i]=condition;
+}
+
+const vector<pair<Plane,Plane> >  Box::getBoundingPlanes(){
+  vector<pair<Plane,Plane> > res;
+  for(int i=0;i<NDIM;++i){
+    Vector base1,base2,base3;
+    base2[(i+1)%NDIM] = 1.;
+    base3[(i+2)%NDIM] = 1.;
+    Plane p1(translateIn(base1),
+             translateIn(base2),
+             translateIn(base3));
+    Vector offset;
+    offset[i]=1;
+    Plane p2(translateIn(base1+offset),
+             translateIn(base2+offset),
+             translateIn(base3+offset));
+    res.push_back(make_pair(p1,p2));
+  }
+  return res;
+}
+
+void Box::setCuboid(const Vector &endpoint){
+  ASSERT(endpoint[0]>0 && endpoint[1]>0 && endpoint[2]>0,"Vector does not define a full cuboid");
+  M->one();
+  M->diagonal()=endpoint;
+  Vector &dinv = Minv->diagonal();
+  for(int i=NDIM;i--;)
+    dinv[i]=1/endpoint[i];
 }
 
@@ -152 +266 @@
     M    = new Matrix(*src.M);
     Minv = new Matrix(*src.Minv);
+    conditions = src.conditions;
   }
   return *this;