Changeset 6c5812 for src/boundary.cpp

Timestamp:

Jun 12, 2008, 7:34:36 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:

318bfd

Parents:

110ceb

Message:

Smaller fixes

CreateClustersinWater is now PrepareClustersinWater (as actual suspension is done in VMD still)
BoundaryPoints can be returned or not in VolumeOfConvexEnvelope

File:

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

src/boundary.cpp

@@ -169 +169 @@
 // ========================================== F U N C T I O N S =================================
 
+/** Finds the endpoint two lines are sharing.
+ * \param *line1 first line
+ * \param *line2 second line
+ * \return point which is shared or NULL if none
+ */
 class BoundaryPointSet * GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
 {
@@ -437 +442 @@
  * \param *out output stream for debugging
  * \param *configuration needed for path to store convex envelope file
+ * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
  * \param *mol molecule structure representing the cluster
  */
-double VolumeOfConvexEnvelope(ofstream *out, config *configuration, molecule *mol) 
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
   atom *Walker = NULL;
   struct Tesselation *TesselStruct = new Tesselation;
-  Boundaries *BoundaryPoints = NULL;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
   
   // 1. calculate center of gravity
@@ -459 +466 @@
   
   // 3. Find all points on the boundary
-  BoundaryPoints = GetBoundaryPoints(out, mol);
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
   
   // 3d. put into boundary set only those points appearing in each of the NDIM sets
@@ -544 +556 @@
 
   // free reference lists
-
+  if (BoundaryFreeFlag)
+    delete[](BoundaryPoints);
+  
   return volume;
 };
@@ -550 +564 @@
 
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
+ * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
  * \param *out output stream for debugging
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
- * \param clustervolume volume of the convex envelope cluster, \sa VolumeOfConvexEnvelope()
+ * \param repetition[] number of repeated cluster per axis
  * \param celldensity desired average density in final cell
- * \param GreatestDiameter[]  greatest diameter of the cluster, \sa GetDiametersOfCluster()
  */
-void CreateClustersinWater(ofstream *out, config *configuration, molecule *mol, double clustervolume, double celldensity, double GreatestDiameter[NDIM])
-{
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity)
+{
+  // 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];
+  int TotalNoClusters = 1;
+  for (int i=0;i<NDIM;i++)
+    TotalNoClusters *= repetition[i];
+
   // sum up the atomic masses
   double totalmass = 0.;
@@ -575 +598 @@
   double cellvolume;
   if (IsAngstroem)
-    cellvolume = (4*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
   else
-    cellvolume = (4*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
+    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 = 4*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
+  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)
     cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
-  double  a = 4*(GreatestDiameter[0]+GreatestDiameter[1]+GreatestDiameter[2]);
-  double b = 4*(GreatestDiameter[0]*GreatestDiameter[1] + GreatestDiameter[0]*GreatestDiameter[2] + GreatestDiameter[1]*GreatestDiameter[2]);
-  double c = minimumvolume - cellvolume;
+  
+  coeffs[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+  coeffs[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;
   double x0 = 0.,x1 = 0.,x2 = 0.;
-  if (gsl_poly_solve_cubic(a,b,c,&x0,&x1,&x2) == 1) // either 1 or 3 on return
+  if (gsl_poly_solve_cubic(coeffs[0],coeffs[1],coeffs[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
     *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
   else {
@@ -593 +619 @@
     x0 = x2;  // sorted in ascending order
   }
-  
-  a = 2 * (x0 + GreatestDiameter[0]);
-  b = 2 * (x0 + GreatestDiameter[1]);
-  c = (x0 + GreatestDiameter[2]);
-  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << a << " and " << b << " and " << c << " with total volume of " << a*b*c << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+
+  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;
 };