Changeset 21c017 for src/molecules.cpp

Timestamp:

Jul 10, 2009, 12:48:05 PM (16 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:

5466f3

Parents:

a37350

Message:

molecule::CenterInBox puts atoms now periodically into the given box, new function molecule::TranslatePeriodically, BUGFIX: molecule::ReturnFullMatrixforSymmetrical()

• molecule::CenterInBox() has no more a vector as a parameter, but instead enforces the periodicity of the simulation box, i.e. all atoms out of bounds are put back in with wrap-around at boundaries. Call of function was changed in everywhere.
• in ParseCommandLineParameters() a SetBoxDimension was missing in certain Center...() commands.
• new function molecule::TranslatePeriodically translates all atoms of a molecule while adhering to the periodicity of the domain
• new function vector::InverseMatrix() returns the hard-encoded inverse of 3x3 real matrix
• BUGFIX: molecule::ReturnFullMatrixforSymmetrical()'s assignment from 6-doubles to 9-doubles was all wrong (symmetric to full 3x3 matrix)

File:

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

src/molecules.cpp

@@ -629 +629 @@
 /** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths.
  * \param *out output stream for debugging
- * \param *BoxLengths box lengths
- */
-bool molecule::CenterInBox(ofstream *out, Vector *BoxLengths)
+ */
+bool molecule::CenterInBox(ofstream *out)
 {
         bool status = true;
         atom *ptr = NULL;
-        Vector *min = new Vector;
-        Vector *max = new Vector;
-
-        // gather min and max for each axis
-        ptr = start->next;      // start at first in list
-        if (ptr != end) {        //list not empty?
-                for (int i=NDIM;i--;) {
-                        max->x[i] = ptr->x.x[i];
-                        min->x[i] = ptr->x.x[i];
-                }
-                while (ptr->next != end) {      // continue with second if present
-                        ptr = ptr->next;
-                        //ptr->Output(1,1,out);
-                        for (int i=NDIM;i--;) {
-                                max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
-                                min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
-                        }
-                }
-        }
-        // sanity check
-        for(int i=NDIM;i--;) {
-                if (max->x[i] - min->x[i] > BoxLengths->x[i])
-                        status = false;
-        }
-        // warn if check failed
-        if (!status)
-                *out << "WARNING: molecule is bigger than defined box!" << endl;
-        else {  // else center in box
-                max->AddVector(min);
-                max->Scale(-1.);
-                max->AddVector(BoxLengths);
-                max->Scale(0.5);
-                Translate(max);
-        }
-
-        // free and exit
-        delete(min);
-        delete(max);
+        Vector x;
+        double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = x.InverseMatrix(M);
+        double value;
+
+//      cout << "The box matrix is :" << endl;
+//  for (int i=0;i<NDIM;++i) {
+//    for (int j=0;j<NDIM;++j)
+//      cout << M[i*NDIM+j] << "\t";
+//    cout << endl;
+//  }
+//  cout << "And its inverse is :" << endl;
+//  for (int i=0;i<NDIM;++i) {
+//    for (int j=0;j<NDIM;++j)
+//      cout << Minv[i*NDIM+j] << "\t";
+//    cout << endl;
+//  }
+        // go through all atoms
+  ptr = start->next;  // start at first in list
+  if (ptr != end) {  //list not empty?
+    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++) {
+        value = floor(x.x[i]);  // next lower integer
+        if (x.x[i] >=0) {
+          x.x[i] -= value;
+        } else {
+          x.x[i] += value+1;
+        }
+      }
+      // matrix multiply
+      x.MatrixMultiplication(M);
+      ptr->x.CopyVector(&x);
+    }
+  }
+  delete(M);
+  delete(Minv);
         return status;
 };
@@ -836 +839 @@
         }
 };
+
+/** Translate the molecule periodically in the box.
+ * \param trans[] translation vector.
+ */
+void molecule::TranslatePeriodically(const Vector *trans)
+{
+  atom *ptr = NULL;
+  Vector x;
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = x.InverseMatrix(M);
+  double value;
+
+  // go through all atoms
+  ptr = start->next;  // start at first in list
+  if (ptr != end) {  //list not empty?
+    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.Translate(trans);
+      x.MatrixMultiplication(Minv);
+      // truncate to [0,1] for each axis
+      for (int i=0;i<NDIM;i++) {
+        value = floor(fabs(x.x[i]));  // next lower integer
+        if (x.x[i] >=0) {
+          x.x[i] -= value;
+        } else {
+          x.x[i] += value+1;
+        }
+      }
+      // matrix multiply
+      x.MatrixMultiplication(M);
+      ptr->x.CopyVector(&x);
+      for (int j=0;j<MDSteps;j++) {
+        x.CopyVector(&Trajectories[ptr].R.at(j));
+        x.Translate(trans);
+        x.MatrixMultiplication(Minv);
+        // truncate to [0,1] for each axis
+        for (int i=0;i<NDIM;i++) {
+          value = floor(x.x[i]);  // next lower integer
+          if (x.x[i] >=0) {
+            x.x[i] -= value;
+          } else {
+            x.x[i] += value+1;
+          }
+        }
+        // matrix multiply
+        x.MatrixMultiplication(M);
+        Trajectories[ptr].R.at(j).CopyVector(&x);
+      }
+    }
+  }
+  delete(M);
+  delete(Minv);
+};
+
 
 /** Mirrors all atoms against a given plane.
@@ -4130 +4190 @@
         matrix[0] = symm[0];
         matrix[1] = symm[1];
-        matrix[2] = symm[3];
+        matrix[2] = symm[2];
         matrix[3] = symm[1];
-        matrix[4] = symm[2];
+        matrix[4] = symm[3];
         matrix[5] = symm[4];
-        matrix[6] = symm[3];
+        matrix[6] = symm[2];
         matrix[7] = symm[4];
         matrix[8] = symm[5];