Changeset f3278b for src

Timestamp:

Aug 7, 2009, 9:13:21 AM (15 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:

1999d8

Parents:

1953f9

Message:

Test case of filling a simulation domain with water included.

Location:

src

Files:

• boundary.cpp (modified) (3 diffs)
• boundary.hpp (modified) (1 diff)
• builder.cpp (modified) (4 diffs)
• molecules.cpp (modified) (1 diff)
• molecules.hpp (modified) (1 diff)

src/boundary.cpp

@@ -863 +863 @@
 
             // insert into Filling
+
+            // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
             *out << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
             Filling->AddAtom(CopyAtoms[Walker->nr]);
@@ -930 +932 @@
       if (!failflag)
         cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+      // write temporary envelope
+      if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % SingleStepWidth == 0))) { // if we have a new triangle and want to output each new triangle configuration
+        TriangleMap::iterator runner = mol->TesselStruct->TrianglesOnBoundary.end();
+        runner--;
+        class BoundaryTriangleSet *triangle = runner->second;
+        if (triangle != NULL) {
+          sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
+          if (DoTecplotOutput) {
+            string NameofTempFile(filename);
+            NameofTempFile.append(NumberName);
+            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+            NameofTempFile.erase(npos, 1);
+            NameofTempFile.append(TecplotSuffix);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
+            tempstream->close();
+            tempstream->flush();
+            delete(tempstream);
+          }
+
+          if (DoRaster3DOutput) {
+            string NameofTempFile(filename);
+            NameofTempFile.append(NumberName);
+            for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+            NameofTempFile.erase(npos, 1);
+            NameofTempFile.append(Raster3DSuffix);
+            *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+            tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+            write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
+    //        // include the current position of the virtual sphere in the temporary raster3d file
+    //        // make the circumsphere's center absolute again
+    //        helper.CopyVector(BaseRay->endpoints[0]->node->node);
+    //        helper.AddVector(BaseRay->endpoints[1]->node->node);
+    //        helper.Scale(0.5);
+    //        (*it)->OptCenter.AddVector(&helper);
+    //        Vector *center = mol->DetermineCenterOfAll(out);
+    //        (*it)->OptCenter.SubtractVector(center);
+    //        delete(center);
+    //        // and add to file plus translucency object
+    //        *tempstream << "# current virtual sphere\n";
+    //        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+    //        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+    //          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+    //          << "\t" << RADIUS << "\t1 0 0\n";
+    //        *tempstream << "9\n  terminating special property\n";
+            tempstream->close();
+            tempstream->flush();
+            delete(tempstream);
+          }
+        }
+        if (DoTecplotOutput || DoRaster3DOutput)
+          TriangleFilesWritten++;
+      }
     } else {
       //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
@@ -941 +997 @@
       baseline = mol->TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
       flag = false;
-    }
-
-    // write temporary envelope
-    if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      class BoundaryTriangleSet *triangle = (mol->TesselStruct->TrianglesOnBoundary.end()--)->second;
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-
-      if (DoRaster3DOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
-//        // include the current position of the virtual sphere in the temporary raster3d file
-//        // make the circumsphere's center absolute again
-//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
-//        helper.AddVector(BaseRay->endpoints[1]->node->node);
-//        helper.Scale(0.5);
-//        (*it)->OptCenter.AddVector(&helper);
-//        Vector *center = mol->DetermineCenterOfAll(out);
-//        (*it)->OptCenter.SubtractVector(center);
-//        delete(center);
-//        // and add to file plus translucency object
-//        *tempstream << "# current virtual sphere\n";
-//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-//          << "\t" << RADIUS << "\t1 0 0\n";
-//        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
     }
   }

src/boundary.hpp

@@ -17 +17 @@
 #define DEBUG 1
 #define DoSingleStepOutput 0
+#define SingleStepWidth 1
 #define DoTecplotOutput 1
 #define DoRaster3DOutput 1

src/builder.cpp

@@ -1348 +1348 @@
             cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
             cout << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
-            cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;
+            cout << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << 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 <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
@@ -1358 +1358 @@
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
+            cout << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl;
             cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
             cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
@@ -1506 +1507 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
-            case 'B':
+            case 'M':
               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                 ExitFlag = 255;
@@ -1767 +1768 @@
                 // center
                 mol->CenterInBox((ofstream *)&cout);
+                argptr+=6;
+              }
+              break;
+            case 'B':
+              ExitFlag = 1;
+              if ((argptr+5 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
+              } else {
+                SaveFlag = true;
+                j = -1;
+                cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+                for (int i=0;i<6;i++) {
+                  mol->cell_size[i] = atof(argv[argptr+i]);
+                }
+                // center
+                mol->BoundInBox((ofstream *)&cout);
                 argptr+=6;
               }

src/molecules.cpp

@@ -741 +741 @@
   delete(Minv);
   delete(Center);
+  return status;
+};
+
+
+/** Bounds the molecule in the box whose lengths are defined by vector \a *BoxLengths.
+ * \param *out output stream for debugging
+ */
+bool molecule::BoundInBox(ofstream *out)
+{
+  bool status = true;
+  Vector x;
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = x.InverseMatrix(M);
+
+  // go through all atoms
+  atom *ptr = start;  // start at first in list
+  while (ptr->next != end) {  // continue with second if present
+    ptr = ptr->next;
+    //ptr->Output(1,1,out);
+    // multiply its vector with matrix inverse
+    x.CopyVector(&ptr->x);
+    x.MatrixMultiplication(Minv);
+    // truncate to [0,1] for each axis
+    for (int i=0;i<NDIM;i++) {
+      while (x.x[i] >= 1.)
+        x.x[i] -= 1.;
+      while (x.x[i] < 0.)
+        x.x[i] += 1.;
+    }
+    x.MatrixMultiplication(M);
+    ptr->x.CopyVector(&x);
+  }
+  delete(M);
+  delete(Minv);
   return status;
 };

src/molecules.hpp

@@ -175 +175 @@
   void CalculateOrbitals(class config &configuration);
   bool CenterInBox(ofstream *out);
+  bool BoundInBox(ofstream *out);
   void CenterEdge(ofstream *out, Vector *max);
   void CenterOrigin(ofstream *out);