Changeset 16648f for src

Timestamp:

Jul 29, 2010, 4:09:57 PM (15 years ago)

Author:

Tillmann Crueger <crueger@…>

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:

04e61a

Parents:

f3be87

Message:

Made the Box::explode() method handle boundary conditions

File:

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

src/Box.cpp

@@ -11 +11 @@
 
 #include <cmath>
+#include <iostream>
 
 #include "Matrix.hpp"
@@ -17 +18 @@
 
 #include "Helpers/Assert.hpp"
+
+using namespace std;
 
 Box::Box()
@@ -105 +108 @@
 
 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;
+  }
+
+  int pos = 0;
+  while(pos<dims){
+    // 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
+    pos=0;
+    ++index[pos];
+    while(index[pos]>n && pos < dims){
+      index[pos++]=-n;
+      ++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);
 }