Changeset 318bfd for src/boundary.cpp

Timestamp:

Jun 13, 2008, 9:18:57 AM (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:

edb650

Parents:

6c5812

Message:

Final touches for Create Clusters in water

repetition in CreateClustersinWater(), some stuff thrown out in VolumeOfConvexEnvelope (loop to throw out some more points, which does not happen if we take any appearing in at least one BoundaryPoints), some variable declarations shifted to the start of a function, BoundaryPoints are given to functions and if NULL is created by calling GetBoundaryPoints() and free'd subsequently.

File:

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

src/boundary.cpp

@@ -372 +372 @@
  * \param *out output stream for debugging
  * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
+ * \param *mol molecule structure representing the cluster
  * \param IsAngstroem whether we have angstroem or atomic units
  * \return NDIM array of the diameters
  */ 
-double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPoints, bool IsAngstroem)
-{
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
+{
+  // get points on boundary of NULL was given as parameter
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+  
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -434 +445 @@
   *out << Verbose(0) << "RESULT: The biggest diameters are " << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "." << endl;
 
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+
   return GreatestDiameter;
 };
@@ -452 +467 @@
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
-  
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G,h;
+  vector x,y;
+  double a,b,c;
+
   // 1. calculate center of gravity
   *out << endl;
@@ -473 +493 @@
   }
   
-  // 3d. put into boundary set only those points appearing in each of the NDIM sets
-  int *AtomList = new int[mol->AtomCount];
-  for (int j=0; j<mol->AtomCount; j++) // reset list
-    AtomList[j] = 0;
-  for (int axis=0; axis<NDIM; axis++)  {  // fill list when it's on the boundary
+  // 4. fill the boundary point list
+  for (int axis=0;axis<NDIM;axis++)
     for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-      AtomList[runner->second.second->nr]++;
+      TesselStruct->AddPoint(runner->second.second);
     }
-  }
-  for (int axis=0; axis<NDIM; axis++)  {
-    for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-      if (AtomList[runner->second.second->nr] < 1) {
-        *out << Verbose(1) << "Throwing especially out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
-        BoundaryPoints[axis].erase(runner);
-      }
-    }
-  }
-  delete[](AtomList);
-  
-  // 4a. fill the boundary point list
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    if (AtomList[Walker->nr] > 0) {
-      TesselStruct->AddPoint(Walker);
-    }
-  }
 
   *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
@@ -512 +510 @@
 //  *out << endl;
   
-  // 4b. guess starting triangle
+  // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
   
-  // 5. go through all lines, that are not yet part of two triangles (only of one so far)
+  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
   TesselStruct->TesselateOnBoundary(out, configuration, mol);
 
@@ -522 +520 @@
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
   *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
-  double volume = 0.;
-  double PyramidVolume = 0.;
-  double G,h;
-  vector x,y;
-  double a,b,c;
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
     x.CopyVector(&runner->second->endpoints[0]->node->x);
@@ -568 +561 @@
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
- * \param repetition[] number of repeated cluster per axis
  * \param celldensity desired average density in final cell
  */
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity)
-{
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity)
+{
+  // transform to PAS
+  mol->PrincipalAxisSystem(out, true);
+  
   // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
   Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
   double clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
-  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, IsAngstroem);
-  double coeffs[NDIM];
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
+  vector BoxLengths;
+  int repetition[NDIM] = {1, 1, 1};
   int TotalNoClusters = 1;
   for (int i=0;i<NDIM;i++)
@@ -602 +598 @@
     cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
   *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
   double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
   *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  if (minimumvolume < cellvolume)
+  if (minimumvolume > cellvolume)
     cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
   
-  coeffs[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
-  coeffs[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
+  BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+  BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
             + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
             + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
-  coeffs[2] = minimumvolume - cellvolume;
+  BoxLengths.x[2] = minimumvolume - cellvolume;
   double x0 = 0.,x1 = 0.,x2 = 0.;
-  if (gsl_poly_solve_cubic(coeffs[0],coeffs[1],coeffs[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
+  if (gsl_poly_solve_cubic(BoxLengths.x[0],BoxLengths.x[1],BoxLengths.x[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
     *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
   else {
@@ -622 +619 @@
   cellvolume = 1;
   for(int i=0;i<NDIM;i++) {
-    coeffs[i] = repetition[0] * (x0 + GreatestDiameter[0]);
-    cellvolume *= coeffs[i];
-  }
-  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << coeffs[0] << " and " << coeffs[1] << " and " << coeffs[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+    BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+    cellvolume *= BoxLengths.x[i];
+  }
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
+  // set new box dimensions
+  *out << Verbose(0) << "Translating to box with these boundaries." << endl;
+  mol->CenterInBox((ofstream *)&cout, &BoxLengths);
+  // update Box of atoms by boundary
+  mol->SetBoxDimension(&BoxLengths);
+  
 };