Changeset 318bfd for src

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.

Location:

src

Files:

• boundary.cpp (modified) (9 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (3 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);
+  
 };
 

src/boundary.hpp

@@ -123 +123 @@
 
 double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPoints, molecule *mol);
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, int repetition[NDIM], double celldensity);
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double celldensity);
 
 

src/builder.cpp

@@ -758 +758 @@
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
@@ -768 +769 @@
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
             cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl; 
             cout << "\t-v/-V\t\tGives version information." << endl;
             cout << "Note: config files must not begin with '-' !" << endl;
@@ -983 +985 @@
               ExitFlag = 1;
               cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-//              tmp = atof(argv[argptr++]);
-//              if (tmp < 1.0) {
-//                cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
-//                tmp = 1.3;
-//              }
               VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
+              break;
+            case 'u':
+              {
+                double tmp;
+                ExitFlag = 1;
+                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
+                tmp = atof(argv[argptr++]);
+                if (tmp < 1.0) {
+                  cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  tmp = 1.3;
+                }
+//                for(int i=0;i<NDIM;i++) {
+//                  repetition[i] = atoi(argv[argptr++]);
+//                  if (repetition[i] < 1)
+//                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
+//                  repetition[i] = 1;
+//                }
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, tmp);
+              }
+              break;
+            case 'd':
+              for (int axis = 1; axis <= NDIM; axis++) {
+                int faktor = atoi(argv[argptr++]);
+                int count;
+                element ** Elements;
+                vector ** Vectors;
+                if (faktor < 1) {
+                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
+                  faktor = 1;
+                }
+                mol->CountAtoms((ofstream *)&cout);  // recount atoms
+                if (mol->AtomCount != 0) {  // if there is more than none
+                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
+                  Elements = new element *[count];
+                  Vectors = new vector *[count];
+                  j = 0;
+                  first = mol->start;
+                  while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
+                    first = first->next;
+                    Elements[j] = first->type;
+                    Vectors[j] = &first->x;
+                    j++;
+                  }
+                  if (count != j)
+                    cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+                  x.Zero();
+                  y.Zero();
+                  y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+                  for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
+                    x.AddVector(&y); // per factor one cell width further
+                    for (int k=count;k--;) { // go through every atom of the original cell
+                      first = new atom(); // create a new body
+                      first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+                      first->x.AddVector(&x);      // translate the coordinates
+                      first->type = Elements[k];  // insert original element
+                      mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+                    }
+                  }
+                  // free memory
+                  delete[](Elements);
+                  delete[](Vectors);
+                  // correct cell size
+                  if (axis < 0) { // if sign was negative, we have to translate everything
+                    x.Zero();
+                    x.AddVector(&y);
+                    x.Scale(-(faktor-1));
+                    mol->Translate(&x);
+                  }
+                  mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor; 
+                }
+              }
               break;
             default:   // no match? Step on