Changeset b12a35 for src/analyzer.cpp


Ignore:
Timestamp:
Oct 30, 2008, 4:17:23 PM (17 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:
05d2b2
Parents:
f731ae
git-author:
Frederik Heber <heber@…> (10/19/08 13:57:23)
git-committer:
Frederik Heber <heber@…> (10/30/08 16:17:23)
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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analyzer.cpp

    rf731ae rb12a35  
    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
Note: See TracChangeset for help on using the changeset viewer.