Ignore:
Timestamp:
Jul 23, 2009, 9:14:18 AM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
3b470f
Parents:
2218dca (diff), 6ce722 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'HessianMatrix'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/analyzer.cpp

    r2218dca r09f4dc  
    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;
    … …  
    3237  ForceMatrix ChiPAS;
    3338  EnergyMatrix Time;
    34   EnergyMatrix EnergyFragments;
    35   EnergyMatrix HcorrectionFragments;
    36   ForceMatrix ForceFragments;
    3739  ForceMatrix ShieldingFragments;
    3840  ForceMatrix ShieldingPASFragments;
    … …  
    5355  stringstream yrange;
    5456  char *dir = NULL;
    55   bool Hcorrected = true;
     57  bool NoHCorrection = false;
     58  bool NoHessian = false;
     59  bool NoTime = false;
    5660  double norm;
    5761  int counter;
    … …  
    8892  // ------------- Parse through all Fragment subdirs --------
    8993  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
    90   Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
     94  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
     95    NoHCorrection = true;
     96    cout << "No HCorrection file found, skipping these." << endl;
     97  }
     98 
    9199  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
    92   if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
     100  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
     101    NoHessian = true;
     102    cout << "No Hessian file found, skipping these." << endl;
     103  }
     104  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
     105    NoTime = true;
     106    cout << "No speed file found, skipping these." << endl;
     107  }
    93108  if (periode != NULL) { // also look for PAS values
    94109    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    … …  
    99114
    100115  // ---------- Parse the TE Factors into an array -----------------
    101   if (!Energy.ParseIndices()) return 1;
    102   if (Hcorrected) Hcorrection.ParseIndices();
     116  if (!Energy.InitialiseIndices()) return 1;
     117  if (!NoHCorrection)
     118    Hcorrection.InitialiseIndices();
    103119 
    104120  // ---------- Parse the Force indices into an array ---------------
    105121  if (!Force.ParseIndices(argv[1])) return 1;
    106122  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    107   if (!ForceFragments.ParseIndices(argv[1])) return 1;
     123  if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     124
     125  // ---------- Parse hessian indices into an array -----------------
     126  if (!NoHessian) {
     127    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     128    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
     129    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     130  }
    108131
    109132  // ---------- Parse the shielding indices into an array ---------------
    … …  
    129152  // ---------- Parse fragment files created by 'joiner' into an array -------------
    130153  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
    131   if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
     154  if (!NoHCorrection)
     155    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
    132156  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
     157  if (!NoHessian)
     158    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
    133159  if (periode != NULL) { // also look for PAS values
    134160    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
    … …  
    144170  filename << argv[3] << "/" << "energy-forces.all";
    145171  output.open(filename.str().c_str(), ios::out);
    146   output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
     172  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
    147173  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    148     for(int k=0;k<Energy.ColumnCounter;k++)
     174    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
    149175      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
    150176    output << endl;
    … …  
    152178  output << endl;
    153179 
    154   output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
     180  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
    155181  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    156     for(int k=0;k<Force.ColumnCounter;k++)
     182    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
    157183      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    158184    output << endl;
    … …  
    160186  output << endl;
    161187
     188  if (!NoHessian) {
     189    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
     190    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     191      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
     192        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     193      output << endl;
     194    }
     195    output << endl;
     196  }
     197
    162198  if (periode != NULL) { // also look for PAS values
    163     output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
     199    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
    164200    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
    165       for(int k=0;k<Shielding.ColumnCounter;k++)
     201      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
    166202        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
    167203      output << endl;
    … …  
    169205    output << endl;
    170206 
    171     output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     207    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
    172208    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    173       for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     209      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
    174210        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
    175211      output << endl;
    … …  
    194230  }
    195231 
    196   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
    197   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    198     for(int k=0;k<Time.ColumnCounter;k++) {
    199       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    200     }
    201     output << endl;
    202   }
    203   output << endl;
     232  if (!NoTime) {
     233    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
     234    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
     235      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     236        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     237      }
     238      output << endl;
     239    }
     240    output << endl;
     241  }
    204242  output.close();
    205   for(int k=0;k<Time.ColumnCounter;k++)
    206     Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
     243  if (!NoTime)
     244    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
     245      Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
    207246
    208247  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
    … …  
    214253  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    215254  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
    216   if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    217   if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
    218   for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    219     for(int k=Time.ColumnCounter;k--;) {
    220       Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
    221     }
    222   counter = 0;
    223   output << "#Order\tFrag.No.\t" << Time.Header << endl;
    224   output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
    225   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    226     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    227       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    228         for(int k=Time.ColumnCounter;k--;) {
    229           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    230         }
    231     counter += KeySet.FragmentsPerOrder[BondOrder];
    232     output << BondOrder+1 << "\t" << counter;
    233     output2 << BondOrder+1 << "\t" << counter;
    234     for(int k=0;k<Time.ColumnCounter;k++) {
    235       output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    236       if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
    237         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
    238       else
    239         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    240     }
    241     output << endl;
    242     output2 << endl;
    243   }
    244   output.close();
    245   output2.close();
     255  if (!NoTime) {
     256    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
     257    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
     258    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     259      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     260        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
     261      }
     262    counter = 0;
     263    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     264    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     265    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     266      for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     267        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     268          for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     269            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     270          }
     271      counter += KeySet.FragmentsPerOrder[BondOrder];
     272      output << BondOrder+1 << "\t" << counter;
     273      output2 << BondOrder+1 << "\t" << counter;
     274      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     275        output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     276        if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
     277          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
     278        else
     279          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     280      }
     281      output << endl;
     282      output2 << endl;
     283    }
     284    output.close();
     285    output2.close();
     286  }
     287
     288  if (!NoHessian) {
     289    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
     290    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;
     291
     292    if (!CreateDataDeltaFrobeniusOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaFrobeniusHessian_xx-Order", "Plot of error between approximated hessian and full hessian in the frobenius norm versus the Bond Order", datum)) return 1;
     293
     294    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
     295    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
     296    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
     297    output << endl << "# Full" << endl;
     298    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     299      output << j << "\t";
     300      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
     301        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     302      output << endl;
     303    }
     304    output.close();
     305  }
    246306
    247307  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
    … …  
    254314    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    255315      output << j << "\t";
    256       for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     316      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
    257317        output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
    258318      output << endl;
    … …  
    297357  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    298358    output << j << "\t";
    299     for(int k=0;k<Force.ColumnCounter;k++)
     359    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
    300360      output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    301361    output << endl;
    … …  
    346406  yrange.str("[1e-8:1e+1]");
    347407 
    348   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    349   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;
     408  if (!NoTime) {
     409    // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     410    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;
     411  }
    350412 
    351413  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
Note: See TracChangeset for help on using the changeset viewer.