Changeset edb650 for src/boundary.cpp

Timestamp:

Jun 13, 2008, 2:26:50 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:

14d4d4

Parents:

318bfd

Message:

some mor changes to fine-tune the cluster generation.

Sadly, VolumeOfConvexEnvelope() is not robust/stable, with Ratio1.4-1/1800K/compundH.*.xyz no usable convex envelope was found but points were connected right through the cluster.

File:

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

src/boundary.cpp

@@ -561 +561 @@
  * \param *configuration needed for path to store convex envelope file
  * \param *mol molecule structure representing the cluster
+ * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
  * \param celldensity desired average density in final cell
  */
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity)
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity)
 {
   // transform to PAS
@@ -571 +572 @@
   bool IsAngstroem = configuration->GetIsAngstroem();
   Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
-  double clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  double clustervolume;
+  if (ClusterVolume == 0)
+    clustervolume = VolumeOfConvexEnvelope(out, configuration, BoundaryPoints, mol);
+  else 
+    clustervolume = ClusterVolume;
   double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
   vector BoxLengths;
@@ -601 +606 @@
   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;
-  
-  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]);
-  BoxLengths.x[2] = minimumvolume - cellvolume;
-  double x0 = 0.,x1 = 0.,x2 = 0.;
-  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 {
-    *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
-    x0 = x2;  // sorted in ascending order
-  }
-
-  cellvolume = 1;
-  for(int i=0;i<NDIM;i++) {
-    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);
+    cout << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
+    for(int i=0;i<NDIM;i++)
+      BoxLengths.x[i] = GreatestDiameter[i];
+    mol->CenterEdge(out, &BoxLengths);
+  } else {
+    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]);
+    BoxLengths.x[2] = minimumvolume - cellvolume;
+    double x0 = 0.,x1 = 0.,x2 = 0.;
+    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 {
+      *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
+      x0 = x2;  // sorted in ascending order
+    }
+  
+    cellvolume = 1;
+    for(int i=0;i<NDIM;i++) {
+      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths.x[i];
+    }
+  
+    // 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);
-  
+  *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;
 };