Changeset 390248 for src/datacreator.cpp


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

  • src/datacreator.cpp

    rc685eb r390248  
    3030
    3131/** Plots an energy vs. order.
    32  * \param Energy EnergyMatrix class containing matrix values (last matrix is created by summing over Fragments)
    33  * \param EnergyFragments EnergyMatrix class containing matrix values
     32 * \param &Fragments EnergyMatrix class containing matrix values
    3433 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
    3534 * \param *prefix prefix in filename (without ending)
    … …  
    3736 * \return true if file was written successfully
    3837 */
    39 bool CreateDataEnergyOrder(class MatrixContainer &Energy, class MatrixContainer &EnergyFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     38bool CreateDataEnergyOrder(class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
    4039{
    4140  stringstream filename;
    … …  
    4645  cout << msg << endl;
    4746  output << "# " << msg << ", created on " << datum;
    48   output << "#Order\tFrag.No.\t" << Energy.Header << endl;
     47  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
    4948  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    5049    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
    51       for(int j=Energy.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
    52         for(int k=Energy.ColumnCounter;k--;)
    53           Energy.Matrix[Energy.MatrixCounter][j][k] -= EnergyFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     50      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
     51        for(int k=Fragments.ColumnCounter;k--;)
     52          Fragments.Matrix[Fragments.MatrixCounter][j][k] += Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    5453    }
    5554    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    56     for (int l=0;l<Energy.ColumnCounter;l++)
    57       if (fabs(EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]) < MYEPSILON)
    58         output << scientific << "\t" << Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l];
     55    for (int l=0;l<Fragments.ColumnCounter;l++)
     56      output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
     57    output << endl;
     58  }
     59  output.close();
     60  return true;
     61};
     62
     63/** Plots an energy error vs. order.
     64 * \param &Energy EnergyMatrix class containing reference values (in MatrixCounter matrix)
     65 * \param &Fragments EnergyMatrix class containing matrix values
     66 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     67 * \param *prefix prefix in filename (without ending)
     68 * \param *msg message to be place in first line as a comment
     69 * \return true if file was written successfully
     70 */
     71bool CreateDataDeltaEnergyOrder(class EnergyMatrix &Energy, class EnergyMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     72{
     73  stringstream filename;
     74  ofstream output;
     75
     76  filename << prefix << ".dat";
     77  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     78  cout << msg << endl;
     79  output << "# " << msg << ", created on " << datum;
     80  output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
     81  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
     82  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     83    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
     84      for(int j=Fragments.RowCounter[ KeySet.OrderSet[BondOrder][i] ];j--;)
     85        for(int k=Fragments.ColumnCounter;k--;)
     86          Fragments.Matrix[Fragments.MatrixCounter][j][k] -= Fragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     87    }
     88    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     89    for (int l=0;l<Fragments.ColumnCounter;l++)
     90      if (fabs(Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]) < MYEPSILON)
     91        output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l];
    5992      else
    60         output << scientific << "\t" << (Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l] / EnergyFragments.Matrix[EnergyFragments.MatrixCounter][ EnergyFragments.RowCounter[EnergyFragments.MatrixCounter]-1 ][l]);
     93        output << scientific << "\t" << (Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter]-1 ][l] / Energy.Matrix[Energy.MatrixCounter][ Energy.RowCounter[Energy.MatrixCounter]-1 ][l]);
    6194    output << endl;
    6295  }
    … …  
    6699
    67100/** Plot forces vs. order.
    68  */
    69 bool CreateDataForcesOrder(class MatrixContainer &Force, class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
     101 * \param &Fragments ForceMatrix class containing matrix values
     102 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     103 * \param *prefix prefix in filename (without ending)
     104 * \param *msg message to be place in first line as a comment
     105 * \param *datum current date and time
     106 * \return true if file was written successfully
     107 */
     108bool CreateDataForcesOrder(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
    70109{
    71110  stringstream filename;
    … …  
    76115  cout << msg << endl;
    77116  output << "# " << msg << ", created on " << datum;
    78   output << "# Order\tFrag.No.\t" << Force.Header << endl;
    79   for(int j=0;j<=Force.RowCounter[ Force.MatrixCounter ];j++)
    80     for(int k=0;k<Force.ColumnCounter;k++)
    81        Force.Matrix[Force.MatrixCounter][j][k] = 0.;
    82   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    83     for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
    84       for(int l=0;l<Force.RowCounter[ KeySet.OrderSet[BondOrder][i] ];l++) {
    85         int j = Force.Indices[ KeySet.OrderSet[BondOrder][i] ][l];
    86         if (j > Force.RowCounter[Force.MatrixCounter]) {
    87           cerr << "Current force index " << j << " is greater than " << Force.RowCounter[Force.MatrixCounter] << "!" << endl;
    88           return 1;
    89         }
    90         if (j != -1)
    91           for(int k=Force.ColumnCounter;k--;) {
    92             Force.Matrix[Force.MatrixCounter][j][k] += ForceFragments.Matrix[ KeySet.OrderSet[BondOrder][i] ][l][k];
    93           }
     117  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
     118  Fragments.SetLastMatrix(0.,0);
     119  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     120    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
     121    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     122    CreateForce(Fragments, Fragments.MatrixCounter);
     123    for (int l=0;l<Fragments.ColumnCounter;l++)
     124       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     125    output << endl;
     126  }
     127  output.close();
     128  return true;
     129};
     130
     131/** Plot forces error vs. order.
     132 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
     133 * \param &Fragments ForceMatrix class containing matrix values
     134 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     135 * \param *prefix prefix in filename (without ending)
     136 * \param *msg message to be place in first line as a comment
     137 * \param *datum current date and time
     138 * \return true if file was written successfully
     139 */
     140bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int))
     141{
     142  stringstream filename;
     143  ofstream output;
     144
     145  filename << prefix << ".dat";
     146  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     147  cout << msg << endl;
     148  output << "# " << msg << ", created on " << datum;
     149  output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
     150  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
     151  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     152    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
     153    output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
     154    CreateForce(Fragments, Fragments.MatrixCounter);
     155    for (int l=0;l<Fragments.ColumnCounter;l++)
     156       output << scientific << "\t" << Fragments.Matrix[Fragments.MatrixCounter][ Fragments.RowCounter[Fragments.MatrixCounter] ][l];
     157    output << endl;
     158  }
     159  output.close();
     160  return true;
     161};
     162
     163/** Plot forces error vs. vs atom vs. order.
     164 * \param &Force ForceMatrix containing reference values (in MatrixCounter matrix)
     165 * \param &Fragments ForceMatrix class containing matrix values
     166 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     167 * \param *prefix prefix in filename (without ending)
     168 * \param *msg message to be place in first line as a comment
     169 * \param *datum current date and time
     170 * \return true if file was written successfully
     171 */
     172bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     173{
     174  stringstream filename;
     175  ofstream output;
     176  double norm = 0.;
     177
     178  filename << prefix << ".dat";
     179  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     180  cout << msg << endl;
     181  output << "# " << msg << ", created on " << datum;
     182  output << "# AtomNo\t" << Fragments.Header << endl;
     183  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
     184  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     185    //cout << "Current order is " << BondOrder << "." << endl;
     186    Fragments.SumSubForces(Fragments, KeySet, BondOrder, -1.);
     187    // errors per atom
     188    output << "#Order\t" << BondOrder+1 << endl;
     189    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     190      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     191      for (int l=0;l<Fragments.ColumnCounter;l++) {
     192        if (((l+1) % 3) == 0) {
     193          norm = 0.;
     194          for (int m=0;m<NDIM;m++)
     195            norm += Force.Matrix[Force.MatrixCounter][ j ][l+m]*Force.Matrix[Force.MatrixCounter][ j ][l+m];
     196          norm = sqrt(norm);
     197        }                                                                                                           
     198//        if (norm < MYEPSILON)
     199          output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     200//        else
     201//          output << scientific << (Fragments.Matrix[Fragments.MatrixCounter][ j ][l] / norm) << "\t";
    94202      }
    95     }
    96     output << BondOrder+1 << "\t" << KeySet.FragmentsPerOrder[BondOrder];
    97     CreateForce(Force, Force.MatrixCounter);
    98     for (int l=0;l<Force.ColumnCounter;l++)
    99        output << scientific << "\t" << Force.Matrix[Force.MatrixCounter][ Force.RowCounter[Force.MatrixCounter] ][l];
     203      output << endl;
     204    }
     205    output << endl;
     206  }
     207  output.close();
     208  return true;
     209};
     210
     211/** Plot forces error vs. vs atom vs. order.
     212 * \param &Fragments ForceMatrix class containing matrix values
     213 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     214 * \param *prefix prefix in filename (without ending)
     215 * \param *msg message to be place in first line as a comment
     216 * \param *datum current date and time
     217 * \return true if file was written successfully
     218 */
     219bool CreateDataForcesOrderPerAtom(class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     220{
     221  stringstream filename;
     222  ofstream output;
     223
     224  filename << prefix << ".dat";
     225  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     226  cout << msg << endl;
     227  output << "# " << msg << ", created on " << datum;
     228  output << "# AtomNo\t" << Fragments.Header << endl;
     229  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     230    //cout << "Current order is " << BondOrder << "." << endl;
     231    Fragments.SumSubForces(Fragments, KeySet, BondOrder, 1.);
     232    // errors per atom
     233    output << "#Order\t" << BondOrder+1 << endl;
     234    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     235      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     236      for (int l=0;l<Fragments.ColumnCounter;l++)
     237        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     238      output << endl;
     239    }
    100240    output << endl;
    101241  }
    … …  
    286426};
    287427
     428/** Leaves the Force.Matrix as it is.
     429 * \param Force ForceMatrix class containing matrix values
     430 * \param MatrixNumber the index for the ForceMatrix::matrix array
     431*/
     432void CreateSameForce(class MatrixContainer &Force, int MatrixNumber)
     433{
     434  // does nothing
     435};
     436
    288437/** Adds vectorwise all forces.
    289438 * Results are stored in the matrix ForceMatrix::MatrixCounter of \a Force.
    … …  
    295444  for (int l=Force.ColumnCounter;l--;)
    296445    Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] = 0.;
    297   for (int l=5;l<Force.ColumnCounter;l++) {
     446  for (int l=0;l<Force.ColumnCounter;l++) {
    298447    for (int k=Force.RowCounter[MatrixNumber];k--;)
    299448      Force.Matrix[MatrixNumber][ Force.RowCounter[MatrixNumber] ][l] += Force.Matrix[MatrixNumber][k][l];
    … …  
    375524
    376525  line >> item;
    377   for (int i=3; i< Energy.ColumnCounter;i++) {
    378     line >> item;
    379     output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+1 << ")) " << uses;
    380     if (i != (Energy.ColumnCounter-1))
     526  for (int i=2; i<= Energy.ColumnCounter;i++) {
     527    line >> item;
     528    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":(abs($" << i+2 << ")) " << uses;
     529    if (i != (Energy.ColumnCounter))
    381530      output << ", \\";
    382531    output << endl;
    … …  
    397546
    398547  line >> item;
    399   for (int i=3; i< Energy.ColumnCounter;i++) {
    400     line >> item;
    401     output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+1 << " " << uses;
    402     if (i != (Energy.ColumnCounter-1))
     548  for (int i=2; i<= Energy.ColumnCounter;i++) {
     549    line >> item;
     550    output << "'" << prefix << ".dat' title '" << item << "' using " << xargument << ":" << i+2 << " " << uses;
     551    if (i != (Energy.ColumnCounter))
    403552      output << ", \\";
    404553    output << endl;
Note: See TracChangeset for help on using the changeset viewer.