Changeset eeec8f for src


Ignore:
Timestamp:
Oct 19, 2008, 1:57:23 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:
84fc91
Parents:
f731ae
Message:

HessianMatrix implemented fully, but not yet working, probably due to wrong matrix generation in script file (convertHessian.py)

HessianMatrix::IsSymmetric was though to be needed, but is NOT. As we regard full matrices, we don't need to add onto mirrored indices as well
Joiner and Analyzer have seen some small changes and bugfixes: NoHessian was not also always looked at when needed and so on

Location:
src
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • TabularUnified src/analyzer.cpp

    rf731ae reeec8f  
    2525  periodentafel *periode = NULL; // and a period table of all elements
    2626  EnergyMatrix Energy;
     27  EnergyMatrix EnergyFragments;
     28  ForceMatrix Force;
     29  ForceMatrix ForceFragments;
     30  HessianMatrix Hessian;
     31  HessianMatrix HessianFragments;
    2732  EnergyMatrix Hcorrection;
    28   ForceMatrix Force;
     33  EnergyMatrix HcorrectionFragments;
    2934  ForceMatrix Shielding;
    3035  ForceMatrix ShieldingPAS;
    3136  EnergyMatrix Time;
    32   EnergyMatrix EnergyFragments;
    33   EnergyMatrix HcorrectionFragments;
    34   ForceMatrix ForceFragments;
    3537  ForceMatrix ShieldingFragments;
    3638  ForceMatrix ShieldingPASFragments;
     
    4951  stringstream yrange;
    5052  char *dir = NULL;
    51   bool Hcorrected = true;
     53  bool NoHCorrection = false;
     54  bool NoHessian = false;
     55  bool NoTime = false;
    5256  double norm;
    5357  int counter;
     
    8488  // ------------- Parse through all Fragment subdirs --------
    8589  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
    86   Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
     90  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
     91    NoHCorrection = true;
     92    cout << "No HCorrection file found, skipping these." << endl;
     93  }
     94 
    8795  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
    88   if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
     96  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
     97    NoHessian = true;
     98    cout << "No Hessian file found, skipping these." << endl;
     99  }
     100  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
     101    NoTime = true;
     102    cout << "No speed file found, skipping these." << endl;
     103  }
    89104  if (periode != NULL) { // also look for PAS values
    90105    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     
    94109  // ---------- Parse the TE Factors into an array -----------------
    95110  if (!Energy.InitialiseIndices()) return 1;
    96   if (Hcorrected) Hcorrection.InitialiseIndices();
     111  if (!NoHCorrection)
     112    Hcorrection.InitialiseIndices();
    97113 
    98114  // ---------- Parse the Force indices into an array ---------------
    99115  if (!Force.ParseIndices(argv[1])) return 1;
    100116  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    101   if (!ForceFragments.ParseIndices(argv[1])) return 1;
     117  if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     118
     119  // ---------- Parse hessian indices into an array -----------------
     120  if (!NoHessian) {
     121    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     122    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
     123    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     124  }
    102125
    103126  // ---------- Parse the shielding indices into an array ---------------
     
    117140  // ---------- Parse fragment files created by 'joiner' into an array -------------
    118141  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
    119   if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
     142  if (!NoHCorrection)
     143    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
    120144  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
     145  if (!NoHessian)
     146    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
    121147  if (periode != NULL) { // also look for PAS values
    122148    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
     
    130156  filename << argv[3] << "/" << "energy-forces.all";
    131157  output.open(filename.str().c_str(), ios::out);
    132   output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
     158  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
    133159  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    134160    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
     
    138164  output << endl;
    139165 
    140   output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
     166  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
    141167  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    142168    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
     
    146172  output << endl;
    147173
    148   if (periode != NULL) { // also look for PAS values
    149     output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
     174  if (!NoHessian) {
     175    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
     176    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     177      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
     178        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     179      output << endl;
     180    }
     181    output << endl;
     182  }
     183
     184  if (periode != NULL) { // also look for PAS values
     185    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
    150186    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
    151187      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
     
    155191    output << endl;
    156192 
    157     output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     193    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
    158194    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    159195      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
     
    164200  }
    165201 
    166   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
    167   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    168     for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
    169       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    170     }
    171     output << endl;
    172   }
    173   output << endl;
     202  if (!NoTime) {
     203    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
     204    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
     205      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     206        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     207      }
     208      output << endl;
     209    }
     210    output << endl;
     211  }
    174212  output.close();
    175   for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
    176     Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
     213  if (!NoTime)
     214    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
     215      Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
    177216
    178217  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
     
    184223  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    185224  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
    186   if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    187   if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
    188   for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    189     for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
    190       Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
    191     }
    192   counter = 0;
    193   output << "#Order\tFrag.No.\t" << Time.Header << endl;
    194   output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
    195   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    196     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    197       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    198         for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
    199           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    200         }
    201     counter += KeySet.FragmentsPerOrder[BondOrder];
    202     output << BondOrder+1 << "\t" << counter;
    203     output2 << BondOrder+1 << "\t" << counter;
    204     for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
    205       output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    206       if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
    207         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
    208       else
    209         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    210     }
    211     output << endl;
    212     output2 << endl;
    213   }
    214   output.close();
    215   output2.close();
     225  if (!NoTime) {
     226    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
     227    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
     228    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     229      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     230        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
     231      }
     232    counter = 0;
     233    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     234    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     235    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     236      for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     237        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     238          for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     239            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     240          }
     241      counter += KeySet.FragmentsPerOrder[BondOrder];
     242      output << BondOrder+1 << "\t" << counter;
     243      output2 << BondOrder+1 << "\t" << counter;
     244      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     245        output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     246        if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
     247          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
     248        else
     249          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     250      }
     251      output << endl;
     252      output2 << endl;
     253    }
     254    output.close();
     255    output2.close();
     256  }
     257
     258  if (!NoHessian) {
     259    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
     260    if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1;
     261
     262    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
     263    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
     264    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
     265    output << endl << "# Full" << endl;
     266    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     267      output << j << "\t";
     268      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
     269        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     270      output << endl;
     271    }
     272    output.close();
     273  }
    216274
    217275  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
     
    305363  yrange.str("[1e-8:1e+1]");
    306364 
    307   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    308   if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
     365  if (!NoTime) {
     366    // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     367    if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "",  1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1;
     368  }
    309369 
    310370  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
  • TabularUnified src/datacreator.cpp

    rf731ae reeec8f  
    6363  cout << msg << endl;
    6464  output << "# " << msg << ", created on " << datum;
    65   output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
     65  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    6666  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    6767    for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) {
     
    9696  cout << msg << endl;
    9797  output << "# " << msg << ", created on " << datum;
    98   output << "#Order\tFrag.No.\t" << Fragments.Header << endl;
     98  output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    9999  Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0);
    100100  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     
    133133  cout << msg << endl;
    134134  output << "# " << msg << ", created on " << datum;
    135   output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
     135  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    136136  Fragments.SetLastMatrix(0.,0);
    137137  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     
    165165  cout << msg << endl;
    166166  output << "# " << msg << ", created on " << datum;
    167   output << "# Order\tFrag.No.\t" << Fragments.Header << endl;
     167  output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    168168  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0);
    169169  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     
    198198  cout << msg << endl;
    199199  output << "# " << msg << ", created on " << datum;
    200   output << "# AtomNo\t" << Fragments.Header << endl;
     200  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    201201  Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0);
    202202  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     
    244244  cout << msg << endl;
    245245  output << "# " << msg << ", created on " << datum;
    246   output << "# AtomNo\t" << Fragments.Header << endl;
     246  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
    247247  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    248248    //cout << "Current order is " << BondOrder << "." << endl;
     
    262262};
    263263
     264
     265/** Plot hessian error vs. vs atom vs. order.
     266 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix)
     267 * \param &Fragments HessianMatrix class containing matrix values
     268 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     269 * \param *prefix prefix in filename (without ending)
     270 * \param *msg message to be place in first line as a comment
     271 * \param *datum current date and time
     272 * \return true if file was written successfully
     273 */
     274bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     275{
     276  stringstream filename;
     277  ofstream output;
     278  double norm = 0.;
     279
     280  filename << prefix << ".dat";
     281  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     282  cout << msg << endl;
     283  output << "# " << msg << ", created on " << datum;
     284  output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl;
     285  Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0);
     286  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     287    //cout << "Current order is " << BondOrder << "." << endl;
     288    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.);
     289    // errors per atom
     290    output << endl << "#Order\t" << BondOrder+1 << endl;
     291    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     292      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     293      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) {
     294        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     295      }
     296      output << endl;
     297    }
     298    output << endl;
     299  }
     300  output.close();
     301  return true;
     302};
     303
     304/** Plot hessian error vs. vs atom vs. order.
     305 * \param &Fragments HessianMatrix class containing matrix values
     306 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order
     307 * \param *prefix prefix in filename (without ending)
     308 * \param *msg message to be place in first line as a comment
     309 * \param *datum current date and time
     310 * \return true if file was written successfully
     311 */
     312bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum)
     313{
     314  stringstream filename;
     315  ofstream output;
     316
     317  filename << prefix << ".dat";
     318  if (!OpenOutputFile(output, dir, filename.str().c_str())) return false;
     319  cout << msg << endl;
     320  output << "# " << msg << ", created on " << datum;
     321  output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl;
     322  Fragments.SetLastMatrix(0., 0);
     323  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     324    //cout << "Current order is " << BondOrder << "." << endl;
     325    Fragments.SumSubHessians(Fragments, KeySet, BondOrder, 1.);
     326    // errors per atom
     327    output << endl << "#Order\t" << BondOrder+1 << endl;
     328    for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) {
     329      output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t";
     330      for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++)
     331        output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t";
     332      output << endl;
     333    }
     334    output << endl;
     335  }
     336  output.close();
     337  return true;
     338};
     339
    264340/** Plot matrix vs. fragment.
    265341 */
     
    273349  cout << msg << endl;
    274350  output << "# " << msg << ", created on " << datum << endl;
    275   output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
     351  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
    276352  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    277353    for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) {
     
    338414  cout << msg << endl;
    339415  output << "# " << msg << ", created on " << datum;
    340   output << "#Order\tFrag.No.\t" << Fragment.Header << endl;
     416  output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl;
    341417  // max
    342418  for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     
    538614void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
    539615{
    540   stringstream line(Energy.Header);
     616  stringstream line(Energy.Header[ Energy.MatrixCounter ]);
    541617  string token;
    542618
     
    562638void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses)
    563639{
    564   stringstream line(Energy.Header);
     640  stringstream line(Energy.Header[Energy.MatrixCounter]);
    565641  string token;
    566642
     
    586662void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    587663{
    588   stringstream line(Force.Header);
     664  stringstream line(Force.Header[Force.MatrixCounter]);
    589665  string token;
    590666
     
    617693void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    618694{
    619   stringstream line(Force.Header);
     695  stringstream line(Force.Header[Force.MatrixCounter]);
    620696  string token;
    621697
     
    648724void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    649725{
    650   stringstream line(Force.Header);
     726  stringstream line(Force.Header[Force.MatrixCounter]);
    651727  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
    652728  string token;
     
    680756void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses)
    681757{
    682   stringstream line(Force.Header);
     758  stringstream line(Force.Header[Force.MatrixCounter]);
    683759  char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"};
    684760  string token;
  • TabularUnified src/datacreator.hpp

    rf731ae reeec8f  
    2626bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
    2727bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum);
     28bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum);
     29bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum);
    2830bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int));
    2931bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int));
  • TabularUnified src/joiner.cpp

    rf731ae reeec8f  
    3737  stringstream prefix;
    3838  char *dir = NULL;
    39   bool Hcorrected = true;
     39  bool NoHCorrection = false;
    4040  bool NoHessian = false;
    4141
     
    6868  // ------------- Parse through all Fragment subdirs --------
    6969  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
    70   Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
     70  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0)) {
     71    NoHCorrection = true;
     72    cout << "No HCorrection matrices found, skipping these." << endl;
     73  }
    7174  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
    7275  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) {
    7376    NoHessian = true;
    74     cout << "No hessian matrices found, skipping those.";
     77    cout << "No hessian matrices found, skipping these." << endl;
    7578  }
    7679  if (periode != NULL) { // also look for PAS values
     
    8184  // ---------- Parse the TE Factors into an array -----------------
    8285  if (!Energy.InitialiseIndices()) return 1;
    83   if (Hcorrected) Hcorrection.InitialiseIndices();
     86  if (!NoHCorrection)
     87    Hcorrection.InitialiseIndices();
    8488 
    8589  // ---------- Parse the Force indices into an array ---------------
     
    8791
    8892  // ---------- Parse the Hessian (=force) indices into an array ---------------
    89   if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     93  if (!NoHessian)
     94    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
    9095
    9196  // ---------- Parse the shielding indices into an array ---------------
     
    100105  if (!KeySet.ParseManyBodyTerms()) return 1;
    101106  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
    102   if (Hcorrected)  HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
     107  if (!NoHCorrection) 
     108    HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
    103109  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    104110  if (!NoHessian)
     
    126132    cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
    127133    if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
    128     if (Hcorrected) {
     134    if (!NoHCorrection) {
    129135      HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
    130136      if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
    131       if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
     137      Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
    132138    } else
    133139      if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
     
    159165    if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1;
    160166    // hessian
    161     if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1;
     167    if (!NoHessian)
     168      if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1;
    162169    // shieldings
    163170    if (periode != NULL) { // also look for PAS values
     
    170177  prefix << dir << EnergyFragmentSuffix;
    171178  if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
    172   if (Hcorrected) {
     179  if (!NoHCorrection) {
    173180    prefix.str(" ");
    174181    prefix << dir << HcorrectionFragmentSuffix;
     
    195202  // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds
    196203  if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1;
    197   if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
     204  if (!NoHCorrection) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
    198205  if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
    199206  if (!NoHessian)
  • TabularUnified src/parser.cpp

    rf731ae reeec8f  
    5656MatrixContainer::MatrixContainer() {
    5757  Indices = NULL;
    58   Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
     58  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
    5959  Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
    6060  RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
    6161  ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
     62  Header[0] = NULL;
    6263  Matrix[0] = NULL;
    6364  RowCounter[0] = -1;
     
    9091  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
    9192 
    92   Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
     93  if (Header != NULL)
     94    for(int i=MatrixCounter+1;i--;)
     95      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
     96  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
    9397  Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    9498  Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
     
    120124        return false;
    121125      Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
    122       for(int j=Matrix->RowCounter[i];j--;)
     126      for(int j=Matrix->RowCounter[i];j--;) {
    123127        Indices[i][j] = Matrix->Indices[i][j];
     128        //cout << Indices[i][j] << "\t";
     129      }
     130      //cout << endl;
    124131    }
    125132  }
     
    153160    return false;
    154161  }
    155  
    156   // skip some initial lines
     162
     163  // parse header
     164  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
    157165  for (int m=skiplines+1;m--;)
    158     input.getline(Header, 1023);
     166    input.getline(Header[MatrixNr], 1023);
    159167 
    160168  // scan header for number of columns
    161   line.str(Header);
     169  line.str(Header[MatrixNr]);
    162170  for(int k=skipcolumns;k--;)
    163     line >> Header;
     171    line >> Header[MatrixNr];
    164172  //cout << line.str() << endl;
    165173  ColumnCounter[MatrixNr]=0;
     
    169177  }
    170178  //cout << line.str() << endl;
    171   cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
     179  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
    172180  if (ColumnCounter[MatrixNr] == 0)
    173181    cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     
    183191    }
    184192  }
    185   cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
     193  //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
    186194  if (RowCounter[MatrixNr] == 0)
    187195    cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
     
    189197  // allocate matrix if it's not zero dimension in one direction
    190198  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    191     Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     199    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
    192200 
    193201    // parse in each entry for this matrix
     
    195203    input.seekg(ios::beg);
    196204    for (int m=skiplines+1;m--;)
    197       input.getline(Header, 1023);    // skip header
    198     line.str(Header);
     205      input.getline(Header[MatrixNr], 1023);    // skip header
     206    line.str(Header[MatrixNr]);
    199207    for(int k=skipcolumns;k--;)  // skip columns in header too
    200208      line >> filename;
    201     strncpy(Header, line.str().c_str(), 1023); 
     209    strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
    202210    for(int j=0;j<RowCounter[MatrixNr];j++) {
    203       Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     211      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
    204212      input.getline(filename, 1023);
    205213      stringstream lines(filename);
     
    212220      }
    213221      //cout << endl;
    214       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
     222      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
    215223      for(int j=ColumnCounter[MatrixNr];j--;)
    216224        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
     
    266274 
    267275  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
     276  Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
    268277  Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    269278  RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
     
    271280  for(int i=MatrixCounter+1;i--;) {
    272281    Matrix[i] = NULL;
     282    Header[i] = NULL;
    273283    RowCounter[i] = -1;
    274284    ColumnCounter[i] = -1;
     
    287297
    288298/** Allocates and resets the memory for a number \a MCounter of matrices.
    289  * \param *GivenHeader Header line
     299 * \param **GivenHeader Header line for each matrix
    290300 * \param MCounter number of matrices
    291301 * \param *RCounter number of rows for each matrix
    292  * \param *CCounter number of columns for each matrices
     302 * \param *CCounter number of columns for each matrix
    293303 * \return Allocation successful
    294304 */
    295 bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter)
    296 {
    297   Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
    298   strncpy(Header, GivenHeader, 1023);
    299 
     305bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
     306{
    300307  MatrixCounter = MCounter;
    301   Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    302   RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    303   ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
     308  Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
     309  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
     310  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
     311  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
    304312  for(int i=MatrixCounter+1;i--;) {
     313    Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
     314    strncpy(Header[i], GivenHeader[i], 1023);
    305315    RowCounter[i] = RCounter[i];
    306316    ColumnCounter[i] = CCounter[i];
    307     Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]"); 
     317    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]"); 
    308318    for(int j=RowCounter[i]+1;j--;) {
    309       Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     319      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
    310320      for(int k=ColumnCounter[i];k--;)
    311321        Matrix[i][j][k] = 0.;
     
    385395};
    386396
    387 /** Sums the energy with each factor and put into last element of \a ***Energies.
     397/** Sums the entries with each factor and put into last element of \a ***Matrix.
    388398 * Sums over "E"-terms to create the "F"-terms
    389399 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     
    462472      return false;
    463473    }
    464     output << Header << endl;
     474    output << Header[i] << endl;
    465475    for(int j=0;j<RowCounter[i];j++) {
    466476      for(int k=0;k<ColumnCounter[i];k++)
     
    491501    return false;
    492502  }
    493   output << Header << endl;
     503  output << Header[MatrixCounter] << endl;
    494504  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    495505    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
     
    765775      int j = Indices[ FragmentNr ][l];
    766776      if (j > RowCounter[MatrixCounter]) {
    767         cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
     777        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    768778        return false;
    769779      }
     
    772782          int k = Indices[ FragmentNr ][m];
    773783          if (k > ColumnCounter[MatrixCounter]) {
    774             cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << "!" << endl;
     784            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    775785            return false;
    776786          }
     
    785795};
    786796
     797/** Constructor for class HessianMatrix.
     798 */
     799HessianMatrix::HessianMatrix() : MatrixContainer()
     800{
     801   IsSymmetric = true;
     802}
     803
     804/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
     805 * Sums over "E"-terms to create the "F"-terms
     806 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     807 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     808 * \param Order bond order
     809 * \return true if summing was successful
     810 */
     811bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
     812{
     813  // go through each order
     814  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     815    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     816    // then go per order through each suborder and pick together all the terms that contain this fragment
     817    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     818      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     819        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     820          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     821          // if the fragment's indices are all in the current fragment
     822          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     823            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     824            //cout << "Current row index is " << k << "/" << m << "." << endl;
     825            if (m != -1) { // if it's not an added hydrogen
     826              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     827                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     828                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     829                  m = l;
     830                  break; 
     831                }
     832              }
     833              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
     834              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     835                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     836                return false;
     837              }
     838             
     839              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
     840                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
     841                //cout << "Current column index is " << l << "/" << n << "." << endl;
     842                if (n != -1) { // if it's not an added hydrogen
     843                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
     844                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
     845                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
     846                      n = p;
     847                      break; 
     848                    }
     849                  }
     850                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
     851                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     852                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     853                    return false;
     854                  }
     855                  if (Order == SubOrder) { // equal order is always copy from Energies
     856                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     857                  } else {
     858                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     859                  }
     860                }
     861              }
     862            }
     863            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     864             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     865          }
     866        } else {
     867          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     868        }
     869      }
     870    }
     871   //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;
     872  }
     873 
     874  return true;
     875};
    787876
    788877/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
  • TabularUnified src/parser.hpp

    rf731ae reeec8f  
    4545    double ***Matrix;
    4646    int **Indices;
    47     char *Header;
     47    char **Header;
    4848    int MatrixCounter;
    4949    int *RowCounter;
     
    5656  bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr);
    5757  virtual bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
    58   bool AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter);
     58  bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter);
    5959  bool ResetMatrix();
    6060  double FindMinValue();
     
    8888// ======================================= CLASS HessianMatrix =============================
    8989
    90 class HessianMatrix : public MatrixContainer {
     90class HessianMatrix : public MatrixContainer { 
    9191  public:
     92    HessianMatrix();
     93    //~HessianMatrix();
    9294    bool ParseIndices(char *name);
     95    bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order);
    9396    bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign);
    9497    bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns);
     98  private:
     99    bool IsSymmetric;
    95100};
    96101
Note: See TracChangeset for help on using the changeset viewer.