Changeset 21c017 for src

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)

Location:

src

Files:

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

src/boundary.cpp

@@ -1010 +1010 @@
       // set new box dimensions
       *out << Verbose(0) << "Translating to box with these boundaries." << endl;
-      mol->CenterInBox((ofstream *) &cout, &BoxLengths);
+      mol->SetBoxDimension(&BoxLengths);
+      mol->CenterInBox((ofstream *) &cout);
     }
   // update Box of atoms by boundary

src/builder.cpp

@@ -308 +308 @@
                         }
                         // center
-                        mol->CenterInBox((ofstream *)&cout, &x);
+                        mol->CenterInBox((ofstream *)&cout);
                         // update Box of atoms by boundary
                         mol->SetBoxDimension(&x);
@@ -1245 +1245 @@
                                                 cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
                                                 cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
-                                                cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << 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-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
@@ -1263 +1263 @@
                                                 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
                                                 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
-                                                cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
+            cout << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl;
+                                                cout << "\t-S <file> Store temperatures from the config file in <file>." << 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;
@@ -1486 +1487 @@
                                                         }
                                                         break;
-                                                case 'T':
+                                                case 'S':
                                                         ExitFlag = 1;
                                                         if ((argptr >= argc) || (argv[argptr][0] == '-')) {
                                                                 ExitFlag = 255;
-                                                                cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
+                                                                cerr << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;
                                                         } else {
                                                                 cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
@@ -1533 +1534 @@
                                                         }
                                                         break;
+            case 'T':
+              ExitFlag = 1;
+              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;
+              } else {
+                ExitFlag = 1;
+                SaveFlag = true;
+                cout << Verbose(1) << "Translating all ions periodically to new origin." << endl;
+                for (int i=NDIM;i--;)
+                  x.x[i] = atof(argv[argptr+i]);
+                mol->TranslatePeriodically((const Vector *)&x);
+                argptr+=3;
+              }
+              break;
                                                 case 's':
                                                         ExitFlag = 1;
@@ -1562 +1578 @@
                                                 case 'b':
                                                         ExitFlag = 1;
-                                                        if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                                                        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 centering in box: -b <length_x> <length_y> <length_z>" << endl;
+                                                                cerr << "Not enough or invalid arguments given for centering 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;
-                                                                j=-1;
-                                                                for (int i=0;i<NDIM;i++) {
-                                                                        j += i+1;
-                                                                        x.x[i] = atof(argv[argptr++]);
-                                                                        mol->cell_size[j] += x.x[i]*2.;
+                                                                for (int i=0;i<6;i++) {
+                                                                        mol->cell_size[i] = atof(argv[argptr+i]);
                                                                 }
                                                                 // center
-                                                                mol->CenterInBox((ofstream *)&cout, &x);
-                                                                // update Box of atoms by boundary
-                                                                mol->SetBoxDimension(&x);
+                                                                mol->CenterInBox((ofstream *)&cout);
+                argptr+=6;
                                                         }
                                                         break;
@@ -1602 +1613 @@
                                                                 }
                                                                 mol->Translate((const Vector *)&x);
+                argptr+=3;
                                                         }
                                                         break;
@@ -1610 +1622 @@
                                                         mol->CenterOrigin((ofstream *)&cout, &x);
                                                         mol->SetBoxDimension(&x);
+              argptr+=0;
                                                         break;
                                                 case 'r':

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];

src/molecules.hpp

@@ -252 +252 @@
         void CountElements();
         void CalculateOrbitals(class config &configuration);
-        bool CenterInBox(ofstream *out, Vector *BoxLengths);
+        bool CenterInBox(ofstream *out);
         void CenterEdge(ofstream *out, Vector *max);
         void CenterOrigin(ofstream *out, Vector *max);
         void CenterGravity(ofstream *out, Vector *max);
         void Translate(const Vector *x);
+        void TranslatePeriodically(const Vector *trans);
         void Mirror(const Vector *x);
         void Align(Vector *n);

src/vector.cpp

@@ -413 +413 @@
         for (int i=NDIM;i--;)
                 x[i] = C.x[i];
+};
+
+/** Calculate the inverse of a 3x3 matrix.
+ * \param *matrix NDIM_NDIM array
+ */
+double * Vector::InverseMatrix(double *A)
+{
+  double *B = (double *) Malloc(sizeof(double)*NDIM*NDIM, "Vector::InverseMatrix: *B");
+  double detA = RDET3(A);
+  double detAReci;
+
+  for (int i=0;i<NDIM*NDIM;++i)
+    B[i] = 0.;
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+  }
+  return B;
 };
 

src/vector.hpp

@@ -40 +40 @@
         void Scale(double factor);
         void MatrixMultiplication(double *M);
+        double * InverseMatrix(double *A);
         void InverseMatrixMultiplication(double *M);
         void KeepPeriodic(ofstream *out, double *matrix);