Changeset 390248 for src/moleculelist.cpp

Timestamp:

Jul 24, 2008, 2:00:19 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:

958457

Parents:

c685eb

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:

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

src/moleculelist.cpp

@@ -133 +133 @@
 };
 
+/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
+ * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
+ * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
+ * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
+ * \param *out output stream for debugging
+ * \param *path path to file
+ */
+bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
+{
+  atom *Walker = NULL;
+  atom *Runner = NULL;
+  double ***FitConstant = NULL, **correction = NULL;
+  int a,b;
+  ofstream output;
+  ifstream input;
+  string line;
+  stringstream zeile;
+  double distance;
+  char ParsedLine[1023];
+  double tmp;
+  char *FragmentNumber = NULL;
+
+  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
+  // 0. parse in fit constant files that should have the same dimension as the final energy files
+  // 0a. find dimension of matrices with constants
+  line = path;
+  line.append("/");
+  line += FRAGMENTPREFIX;
+  line += "1";
+  line += FITCONSTANTSUFFIX;
+  input.open(line.c_str());
+  if (input == NULL) {
+    cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+    return false;
+  }
+  a=0;
+  b=-1; // we overcount by one 
+  while (!input.eof()) {
+    input.getline(ParsedLine, 1023);
+    zeile.str(ParsedLine);
+    int i=0;
+    while (!zeile.eof()) {
+      zeile >> distance;
+      i++;
+    }
+    if (i > a) 
+      a = i; 
+    b++;
+  }
+  cout << "I recognized " << a << " columns and " << b << " rows, ";
+  input.close();
+  
+  // 0b. allocate memory for constants
+  FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  // 0c. parse in constants
+  for (int i=0;i<3;i++) {
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    sprintf(ParsedLine, "%d", i+1);
+    line += ParsedLine;
+    line += FITCONSTANTSUFFIX;
+    input.open(line.c_str());
+    if (input == NULL) {
+      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+      return false;
+    }
+    int k = 0,l;
+    while ((!input.eof()) && (k < b)) {
+      input.getline(ParsedLine, 1023);
+      //cout << "Current Line: " << ParsedLine << endl;
+      zeile.str(ParsedLine);
+      zeile.clear();
+      l = 0;
+      while ((!zeile.eof()) && (l < a)) {
+        zeile >> FitConstant[i][l][k];
+        //cout << FitConstant[i][l][k] << "\t";
+        l++;
+      }
+      //cout << endl;
+      k++;
+    }
+    input.close();
+  }
+  for(int k=0;k<3;k++) {
+    cout << "Constants " << k << ":" << endl;
+    for (int j=0;j<b;j++) {
+      for (int i=0;i<a;i++) {
+        cout << FitConstant[k][i][j] << "\t";
+      }
+      cout << endl;
+    }
+    cout << endl;
+  }
+  
+  // 0d. allocate final correction matrix
+  correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  for (int i=a;i--;)
+    correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+        
+  // 1a. go through every molecule in the list
+  for(int i=NumberOfMolecules;i--;) {
+    // 1b. zero final correction matrix
+    for (int k=a;k--;)
+      for (int j=b;j--;)
+        correction[k][j] = 0.;
+    // 2. take every hydrogen that is a saturated one
+    Walker = ListOfMolecules[i]->start;
+    while (Walker->next != ListOfMolecules[i]->end) {
+      Walker = Walker->next;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
+        Runner = ListOfMolecules[i]->start;
+        while (Runner->next != ListOfMolecules[i]->end) {
+          Runner = Runner->next;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {  // (hydrogens have only one bonding partner!)
+            // 4. evaluate the morse potential for each matrix component and add up
+            distance = sqrt(Runner->x.Distance(&Walker->x));
+            //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            for(int k=0;k<a;k++) {
+              for (int j=0;j<b;j++) {
+                switch(k) {
+                case 1:
+                case 7:
+                case 11:
+                  tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
+                  break;
+                default:
+                  tmp = FitConstant[0][k][j] * pow( distance,  FitConstant[1][k][j]) + FitConstant[2][k][j];
+                };
+                correction[k][j] -= tmp;    // ground state is actually lower (disturbed by additional interaction)
+                //cout << tmp << "\t";
+              }
+              //cout << endl;
+            }
+            //cout << endl;
+          }
+        }
+      }
+    }
+    // 5. write final matrix to file
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
+    line += FragmentNumber;
+    delete(FragmentNumber);
+    line += HCORRECTIONSUFFIX;
+    output.open(line.c_str());
+    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+    for (int j=0;j<b;j++) {
+      for(int i=0;i<a;i++)
+        output << correction[i][j] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+  line = path;
+  line.append("/");
+  line += HCORRECTIONSUFFIX;
+  output.open(line.c_str());
+  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+  for (int j=0;j<b;j++) {
+    for(int i=0;i<a;i++)
+      output << 0 << "\t";
+    output << endl;
+  }
+  output.close();
+  // 6. free memory of parsed matrices
+  FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k=0;k<3;k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i=a;i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  cout << "done." << endl;
+  return true;
+};
 
 /** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
@@ -224 +411 @@
     outputFragment.open(FragmentName, ios::out);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
-    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
+    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
     else
@@ -263 +450 @@
     configuration->SetDefaultPath(FragmentName);
     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
-    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
+    if ((intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
       *out << " done." << endl;
     else