Changeset db3ea3 for molecuilder/src/parser.cpp

Timestamp:

Jul 24, 2008, 2:00:19 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

747b10

Parents:

b3eb8e

Message:

BIG change: Hcorrection included and bugfix of analyzer (wrt forces)

• Hcorrection is supplied via molecuilder which gives correction values if hydrogen atoms are to close and thus already feel each other's influence. This correction energy matrix is also parsed ans subtracted from the approximated energy matrix generated from the fragments. This impacts both on joiner, analyzer and molecuilder
• Forces were not working which had multiple causes:
  • first, we stepped back down to absolute errors which rather gives a feel about convergence (1e-6 is simply small and neglegible)
  • there may have been bugs where (either in Force or in ForceFragments) adding up took place, this is cleaned up
  • more routines have been abstracted from analyzer into datacreator functions
    • the order of ions had changed, re-run with current molecuilder, calculation and subsequent joining/analyzing brought correct results finally
• plots were bad still due to wrong axes (two more functions finding max/min values), use of wrong columns in plot file

File:

• molecuilder/src/parser.cpp (modified) (8 diffs)

molecuilder/src/parser.cpp

@@ -146 +146 @@
     file.str(" ");
     if (i != MatrixCounter)  
-      file << name << FRAGMENTPREFIX << FragmentNumber << "/" << prefix << suffix;
+      file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
     else
       file << name << prefix << suffix;
@@ -265 +265 @@
 };
 
+/** Scans all elements of MatrixContainer::Matrix for greatest absolute value.
+ * \return greatest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMaxValue()
+{
+  double max = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) > max)
+          max = fabs(Matrix[i][j][k]);
+  if (fabs(max) < MYEPSILON)
+    max += MYEPSILON;
+ return max;
+};
+
+/** Scans all elements of MatrixContainer::Matrix for smallest absolute value.
+ * \return smallest value of MatrixContainer::Matrix
+ */
+double MatrixContainer::FindMinValue()
+{
+  double min = Matrix[0][0][0];
+  for(int i=MatrixCounter+1;i--;)
+    for(int j=RowCounter[i]+1;j--;)
+      for(int k=ColumnCounter;k--;)
+        if (fabs(Matrix[i][j][k]) < min)
+          min = fabs(Matrix[i][j][k]);
+  if (fabs(min) < MYEPSILON)
+    min += MYEPSILON;
+  return min;
+};
+
 /** Sets all values in the last of MatrixContainer::Matrix to \a value.
  * \param value reset value
@@ -333 +365 @@
               }
             }
-          //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
-            //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
+            //if ((ColumnCounter>1) && (RowCounter[0]-1 >= 1))
+             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
           }
         } else {
@@ -341 +373 @@
       }
     }
-    //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
+   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
   }
   
@@ -428 +460 @@
  * Sums over the "F"-terms in ANOVA decomposition.
  * \param Matrix MatrixContainer with matrices (LevelCounter by ColumnCounter) with all the energies.
+ * \param CorrectionFragments MatrixContainer with hydrogen saturation correction per fragments
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
+ * \parsm sign +1 or -1
  * \return true if summing was successful
  */
-bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
+bool EnergyMatrix::SumSubEnergy(class EnergyMatrix &Fragments, class EnergyMatrix *CorrectionFragments, class KeySetsContainer &KeySet, int Order, double sign)
 {
   // sum energy
-  for(int i=KeySet.FragmentsPerOrder[Order];i--;)
-    for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
-      for(int k=ColumnCounter;k--;)
-        Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
+  if (CorrectionFragments == NULL)
+    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter;k--;)
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k];
+  else
+    for(int i=KeySet.FragmentsPerOrder[Order];i--;)
+      for(int j=RowCounter[ KeySet.OrderSet[Order][i] ];j--;)
+        for(int k=ColumnCounter;k--;)
+          Matrix[MatrixCounter][j][k] += sign*(Fragments.Matrix[ KeySet.OrderSet[Order][i] ][j][k] + CorrectionFragments->Matrix[ KeySet.OrderSet[Order][i] ][j][k]);
   return true;
 };
@@ -492 +532 @@
  * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
  * \param Order bond order
+ *  \param sign +1 or -1
  * \return true if summing was successful
  */
-bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order)
-{
+bool ForceMatrix::SumSubForces(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign)
+{
+  int FragmentNr;
   // sum forces
-  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++)
-    for(int l=0;l<RowCounter[ KeySet.OrderSet[Order][i] ];l++) {
-      int j = Indices[ KeySet.OrderSet[Order][i] ][l];
+  for(int i=0;i<KeySet.FragmentsPerOrder[Order];i++) {
+    FragmentNr = KeySet.OrderSet[Order][i];
+    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
+      int j = Indices[ FragmentNr ][l];
       if (j > RowCounter[MatrixCounter]) {
         cerr << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
         return false;
       }
-      if (j != -1)
+      if (j != -1) {
+        //if (j == 0) cout << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
         for(int k=2;k<ColumnCounter;k++)
-          Matrix[MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[Order][i] ][l][k];
-    }
+          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
+      }
+    }
+  }
   return true;
 };
@@ -571 +617 @@
       KeySets[i][j] = -1;
     FragmentNumber = FixedDigitNumber(FragmentCounter, i);
-    cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
+    //cout << FRAGMENTPREFIX << FragmentNumber << "[" << AtomCounter[i] << "]:";
     Free((void **)&FragmentNumber, "KeySetsContainer::ParseKeySets: *FragmentNumber");
     input.getline(filename, 1023);
@@ -577 +623 @@
     for(int j=0;(j<AtomCounter[i]) && (!line.eof());j++) {
       line >> KeySets[i][j];
-      cout << " " << KeySets[i][j];
-    }
-    cout << endl;
+      //cout << " " << KeySets[i][j];
+    }
+    //cout << endl;
   }
   input.close();