Changeset 36ec71 for src/analyzer.cpp


Ignore:
Timestamp:
Jul 24, 2009, 10:38:32 AM (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, Candidate_v1.7.0, 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:
d30402
Parents:
042f82 (diff), 51c910 (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.
git-author:
Frederik Heber <heber@…> (07/23/09 14:23:32)
git-committer:
Frederik Heber <heber@…> (07/24/09 10:38:32)
Message:

Merge branch 'master' into ConcaveHull

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/bond.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/defs.hpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/joiner.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp

merges:

compilation fixes:

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/analyzer.cpp

    r042f82 r36ec71  
    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;
     36  ForceMatrix Chi;
     37  ForceMatrix ChiPAS;
    3138  EnergyMatrix Time;
    32   EnergyMatrix EnergyFragments;
    33   EnergyMatrix HcorrectionFragments;
    34   ForceMatrix ForceFragments;
    3539  ForceMatrix ShieldingFragments;
    3640  ForceMatrix ShieldingPASFragments;
     41  ForceMatrix ChiFragments;
     42  ForceMatrix ChiPASFragments;
    3743  KeySetsContainer KeySet;
    3844  ofstream output;
    … …  
    4955  stringstream yrange;
    5056  char *dir = NULL;
    51   bool Hcorrected = true;
     57  bool NoHessian = false;
     58  bool NoTime = false;
     59  bool NoHCorrection = true;
    5260  int counter;
    5361
    … …  
    8391  // ------------- Parse through all Fragment subdirs --------
    8492  if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1;
    85   Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0);
     93  if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) {
     94    NoHCorrection = true;
     95    cout << "No HCorrection file found, skipping these." << endl;
     96  }
     97 
    8698  if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1;
    87   if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1;
     99  if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) {
     100    NoHessian = true;
     101    cout << "No Hessian file found, skipping these." << endl;
     102  }
     103  if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) {
     104    NoTime = true;
     105    cout << "No speed file found, skipping these." << endl;
     106  }
    88107  if (periode != NULL) { // also look for PAS values
    89108    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    90109    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
     110    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
     111    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
    91112  }
    92113
    93114  // ---------- Parse the TE Factors into an array -----------------
    94115  if (!Energy.ParseIndices()) return 1;
    95   if (Hcorrected) Hcorrection.ParseIndices();
     116  if (!NoHCorrection) Hcorrection.ParseIndices();
    96117
    97118  // ---------- Parse the Force indices into an array ---------------
    98119  if (!Force.ParseIndices(argv[1])) return 1;
    99120  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    100   if (!ForceFragments.ParseIndices(argv[1])) return 1;
     121  if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     122
     123  // ---------- Parse hessian indices into an array -----------------
     124  if (!NoHessian) {
     125    if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     126    if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1;
     127    if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1;
     128  }
    101129
    102130  // ---------- Parse the shielding indices into an array ---------------
    … …  
    108136    if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
    109137    if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
     138    if(!Chi.ParseIndices(argv[1])) return 1;
     139    if(!ChiPAS.ParseIndices(argv[1])) return 1;
     140    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
     141    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
     142    if(!ChiFragments.ParseIndices(argv[1])) return 1;
     143    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
    110144  }
    111145
    … …  
    116150  // ---------- Parse fragment files created by 'joiner' into an array -------------
    117151  if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1;
    118   if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
     152  if (!NoHCorrection)
     153    HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0);
    119154  if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1;
     155  if (!NoHessian)
     156    if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1;
    120157  if (periode != NULL) { // also look for PAS values
    121158    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
    122159    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
     160    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
     161    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
    123162  }
    124163
    … …  
    129168  filename << argv[3] << "/" << "energy-forces.all";
    130169  output.open(filename.str().c_str(), ios::out);
    131   output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;
     170  output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl;
    132171  for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) {
    133     for(int k=0;k<Energy.ColumnCounter;k++)
     172    for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++)
    134173      output << scientific << Energy.Matrix[ Energy.MatrixCounter ][j][k] << "\t";
    135174    output << endl;
    136175  }
    137176  output << endl;
    138 
    139   output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;
     177 
     178  output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl;
    140179  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    141     for(int k=0;k<Force.ColumnCounter;k++)
     180    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
    142181      output << scientific << Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    143182    output << endl;
    … …  
    145184  output << endl;
    146185
     186  if (!NoHessian) {
     187    output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl;
     188    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     189      for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++)
     190        output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     191      output << endl;
     192    }
     193    output << endl;
     194  }
     195
    147196  if (periode != NULL) { // also look for PAS values
    148     output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl;
     197    output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl;
    149198    for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) {
    150       for(int k=0;k<Shielding.ColumnCounter;k++)
     199      for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++)
    151200        output << scientific << Shielding.Matrix[ Shielding.MatrixCounter ][j][k] << "\t";
    152201      output << endl;
    153202    }
    154203    output << endl;
    155 
    156     output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;
     204 
     205    output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl;
    157206    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    158       for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     207      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
    159208        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
    160209      output << endl;
    161210    }
    162211    output << endl;
    163   }
    164 
    165   output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl;
    166   for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
    167     for(int k=0;k<Time.ColumnCounter;k++) {
    168       output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
    169     }
    170     output << endl;
    171   }
    172   output << endl;
     212
     213    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
     214    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
     215      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
     216        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
     217      output << endl;
     218    }
     219    output << endl;
     220 
     221    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
     222    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
     223      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
     224        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
     225      output << endl;
     226    }
     227    output << endl;
     228  }
     229 
     230  if (!NoTime) {
     231    output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl;
     232    for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) {
     233      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     234        output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t";
     235      }
     236      output << endl;
     237    }
     238    output << endl;
     239  }
    173240  output.close();
    174   for(int k=0;k<Time.ColumnCounter;k++)
    175     Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
     241  if (!NoTime)
     242    for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++)
     243      Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k];
    176244
    177245  // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++
    … …  
    183251  // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    184252  // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order
    185   if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
    186   if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
    187   for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    188     for(int k=Time.ColumnCounter;k--;) {
    189       Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
    190     }
    191   counter = 0;
    192   output << "#Order\tFrag.No.\t" << Time.Header << endl;
    193   output2 << "#Order\tFrag.No.\t" << Time.Header << endl;
    194   for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
    195     for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
    196       for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
    197         for(int k=Time.ColumnCounter;k--;) {
    198           Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
    199         }
    200     counter += KeySet.FragmentsPerOrder[BondOrder];
    201     output << BondOrder+1 << "\t" << counter;
    202     output2 << BondOrder+1 << "\t" << counter;
    203     for(int k=0;k<Time.ColumnCounter;k++) {
    204       output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    205       if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
    206         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
    207       else
    208         output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
    209     }
    210     output << endl;
    211     output2 << endl;
    212   }
    213   output.close();
    214   output2.close();
     253  if (!NoTime) {
     254    if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false;
     255    if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false;
     256    for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     257      for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     258        Time.Matrix[ Time.MatrixCounter ][j][k] = 0.;
     259      }
     260    counter = 0;
     261    output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     262    output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl;
     263    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     264      for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;)
     265        for(int j=Time.RowCounter[Time.MatrixCounter];j--;)
     266          for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) {
     267            Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k];
     268          }
     269      counter += KeySet.FragmentsPerOrder[BondOrder];
     270      output << BondOrder+1 << "\t" << counter;
     271      output2 << BondOrder+1 << "\t" << counter;
     272      for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) {
     273        output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     274        if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON)
     275          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k];
     276        else
     277          output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k];
     278      }
     279      output << endl;
     280      output2 << endl;
     281    }
     282    output.close();
     283    output2.close();
     284  }
     285
     286  if (!NoHessian) {
     287    // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM
     288    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;
     289
     290    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;
     291
     292    // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order
     293    if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1;
     294    if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false;
     295    output << endl << "# Full" << endl;
     296    for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) {
     297      output << j << "\t";
     298      for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++)
     299        output << scientific <<  Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t";
     300      output << endl;
     301    }
     302    output.close();
     303  }
    215304
    216305  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
    217 
    218306  if (periode != NULL) { // also look for PAS values
    219307    if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
    … …  
    223311    for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) {
    224312      output << j << "\t";
    225       for(int k=0;k<ShieldingPAS.ColumnCounter;k++)
     313      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
    226314        output << scientific <<  ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
    227315      output << endl;
    228316    }
    229   }
    230   output.close();
     317    output.close();
     318    if (!CreateDataDeltaForcesOrderPerAtom(ChiPAS, ChiPASFragments, KeySet, argv[3], "DeltaChisPAS-Order", "Plot of error between approximated Chis and full Chis versus the Bond Order", datum)) return 1;
     319    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
     320    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
     321    output << endl << "# Full" << endl;
     322    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
     323      output << j << "\t";
     324      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
     325        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
     326      output << endl;
     327    }
     328    output.close();
     329  }
    231330
    232331
    … …  
    255354  for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) {
    256355    output << j << "\t";
    257     for(int k=0;k<Force.ColumnCounter;k++)
     356    for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++)
    258357      output << scientific <<  Force.Matrix[ Force.MatrixCounter ][j][k] << "\t";
    259358    output << endl;
    … …  
    303402  Fragmentxrange << "[0:" << KeySet.FragmentCounter+1 << "]";
    304403  yrange.str("[1e-8:1e+1]");
    305 
    306   // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
    307   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;
    308 
     404 
     405  if (!NoTime) {
     406    // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order
     407    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  }
     409 
    309410  // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM
    310411  if (!CreatePlotOrder(Energy, KeySet, argv[3], "DeltaEnergies-Order", 1, "outside", "y", "",  1, 1, "bond order k", "absolute error in energy [Ht]", Orderxrange.str().c_str(), yrange.str().c_str(), "1" , "with linespoints", AbsEnergyPlotLine)) return 1;
    … …  
    403504    }
    404505    output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    405     output.close();
    406     output2.close();
     506    output2.close(); 
     507
     508    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
     509    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
     510    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     511    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     512    output << "set boxwidth " << step << endl;
     513    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     514    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     515    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     516      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
     517      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
     518      if (BondOrder-1 != KeySet.Order)
     519        output2 << ", \\" << endl;
     520    }
     521    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
     522    output.close(); 
     523    output2.close(); 
     524
     525    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
     526    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
     527    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     528    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     529    output << "set boxwidth " << step << endl;
     530    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     531    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     532    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     533      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
     534      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
     535      if (BondOrder-1 != KeySet.Order)
     536        output2 << ", \\" << endl;
     537    }
     538    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
     539    output.close(); 
     540    output2.close(); 
    407541  }
    408542
Note: See TracChangeset for help on using the changeset viewer.