Changeset 348b6a


Ignore:
Timestamp:
Oct 30, 2008, 12:55:19 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
25f5089
Parents:
5fd79e
Message:

Implemented analysis of magnetic susceptibilites.

This is very similar to Sigma(PAS), however not with ForceMatrix class but with EnergyMatrix. This is still not finished and yet not working. However, it does not impact on other functions of joiner or analyzer.

Location:
molecuilder/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/analyzer.cpp

    r5fd79e r348b6a  
    3434  ForceMatrix Shielding;
    3535  ForceMatrix ShieldingPAS;
     36  ForceMatrix Chi;
     37  ForceMatrix ChiPAS;
    3638  EnergyMatrix Time;
    3739  ForceMatrix ShieldingFragments;
    3840  ForceMatrix ShieldingPASFragments;
     41  ForceMatrix ChiFragments;
     42  ForceMatrix ChiPASFragments;
    3943  KeySetsContainer KeySet;
    4044  ofstream output;
    … …  
    105109    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    106110    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
     111    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
     112    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
    107113  }
    108114
    … …  
    132138    if(!ShieldingFragments.ParseIndices(argv[1])) return 1;
    133139    if(!ShieldingPASFragments.ParseIndices(argv[1])) return 1;
     140    if(!Chi.ParseIndices(argv[1])) return 1;
     141    if(!ChiPAS.ParseIndices(argv[1])) return 1;
     142    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
     143    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
     144    if(!ChiFragments.ParseIndices(argv[1])) return 1;
     145    if(!ChiPASFragments.ParseIndices(argv[1])) return 1;
    134146  }
    135147
    … …  
    148160    if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1;
    149161    if (!ShieldingPASFragments.ParseFragmentMatrix(argv[1], dir, ShieldingPASFragmentSuffix, 1, 0)) return 1;
     162    if (!ChiFragments.ParseFragmentMatrix(argv[1], dir, ChiFragmentSuffix, 1, 0)) return 1;
     163    if (!ChiPASFragments.ParseFragmentMatrix(argv[1], dir, ChiPASFragmentSuffix, 1, 0)) return 1;
    150164  }
    151165
    … …  
    195209      for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++)
    196210        output << scientific << ShieldingPAS.Matrix[ ShieldingPAS.MatrixCounter ][j][k] << "\t";
     211      output << endl;
     212    }
     213    output << endl;
     214
     215    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Hessian.MatrixCounter] << endl;
     216    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
     217      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
     218        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
     219      output << endl;
     220    }
     221    output << endl;
     222 
     223    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
     224    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
     225      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
     226        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
    197227      output << endl;
    198228    }
    … …  
    286316      output << endl;
    287317    }
    288   }
    289   output.close();
     318    output.close();
     319    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;
     320    if (!CreateDataForcesOrderPerAtom(ChiPASFragments, KeySet, argv[3], "ChisPAS-Order", "Plot of approximated Chis versus the Bond Order", datum)) return 1;
     321    if (!AppendOutputFile(output, argv[3], "ChisPAS-Order.dat" )) return false;
     322    output << endl << "# Full" << endl;
     323    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
     324      output << j << "\t";
     325      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
     326        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
     327      output << endl;
     328    }
     329    output.close();
     330  }
    290331
    291332 
    … …  
    466507    output.close(); 
    467508    output2.close(); 
     509
     510    if(!OpenOutputFile(output, argv[3], "ChisPAS-Order.pyx")) return 1;
     511    if(!OpenOutputFile(output2, argv[3], "DeltaChisPAS-Order.pyx")) return 1;
     512    CreatePlotHeader(output, "ChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     513    CreatePlotHeader(output2, "DeltaChisPAS-Order", 1, "top right", NULL, NULL,  1, 5, "nuclei index", "iso chemical Chi value [ppm]");
     514    output << "set boxwidth " << step << endl;
     515    output << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     516    output2 << "plot [0:" << ChiPAS.RowCounter[ChiPAS.MatrixCounter]+10 << "]\\" << endl;
     517    for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) {
     518      output << "'ChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1+" << step*(double)BondOrder << "):7 with boxes, \\" << endl;
     519      output2 << "'DeltaChisPAS-Order.dat' index " << BondOrder << " title 'Order " << BondOrder+1 << "' using ($1):7 with linespoints";
     520      if (BondOrder-1 != KeySet.Order)
     521        output2 << ", \\" << endl;
     522    }
     523    output << "'ChisPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
     524    output.close(); 
     525    output2.close(); 
    468526  }
    469527

  • molecuilder/src/joiner.cpp

    r5fd79e r348b6a  
    3434  ForceMatrix ShieldingFragments;
    3535  ForceMatrix ShieldingPASFragments;
     36  EnergyMatrix Chi;
     37  EnergyMatrix ChiPAS;
     38  EnergyMatrix ChiFragments;
     39  EnergyMatrix ChiPASFragments;
    3640  KeySetsContainer KeySet; 
    3741  stringstream prefix;
    … …  
    8084    if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
    8185    if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
     86    if (!Chi.ParseFragmentMatrix(argv[1], dir, ChiSuffix, 1, 0)) return 1;
     87    if (!ChiPAS.ParseFragmentMatrix(argv[1], dir, ChiPASSuffix, 1, 0)) return 1;
    8288  }
    8389
    … …  
    98104    if(!Shielding.ParseIndices(argv[1])) return 1;
    99105    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
     106    if(!Chi.ParseIndices(argv[1])) return 1;
     107    if(!ChiPAS.ParseIndices(argv[1])) return 1;
    100108  }
    101109
    102110  // ---------- Parse the KeySets into an array ---------------
    103111  if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1;
    104 
    105112  if (!KeySet.ParseManyBodyTerms()) return 1;
     113
    106114  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
    107115  if (!NoHCorrection) 
    … …  
    113121    if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1;
    114122    if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1;
     123    if (!ChiFragments.AllocateMatrix(Chi.Header, Chi.MatrixCounter, Chi.RowCounter, Chi.ColumnCounter)) return 1;
     124    if (!ChiPASFragments.AllocateMatrix(ChiPAS.Header, ChiPAS.MatrixCounter, ChiPAS.RowCounter, ChiPAS.ColumnCounter)) return 1;
    115125  }
    116126 
    … …  
    123133    if(!Shielding.SetLastMatrix(0., 2)) return 1;
    124134    if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
     135    if(!Chi.SetLastMatrix(0., 2)) return 1;
     136    if(!ChiPAS.SetLastMatrix(0., 2)) return 1;
    125137  }
    126138
    … …  
    149161    }
    150162    if (periode != NULL) { // also look for PAS values
    151       cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
     163      cout << "Summing shieldings and susceptibilities of order " << BondOrder+1 << " ..." << endl;
    152164      if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1;
    153165      if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1;
    154166      if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1;
    155167      if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1;
     168      if (!ChiFragments.SumSubManyBodyTerms(Chi, KeySet, BondOrder)) return 1;
     169      if (!Chi.SumSubForces(ChiFragments, KeySet, BondOrder, 1.)) return 1;
     170      if (!ChiPASFragments.SumSubManyBodyTerms(ChiPAS, KeySet, BondOrder)) return 1;
     171      if (!ChiPAS.SumSubForces(ChiPASFragments, KeySet, BondOrder, 1.)) return 1;
    156172    }
    157173
    … …  
    171187      if (!Shielding.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingSuffix)) return 1;
    172188      if (!ShieldingPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingPASSuffix)) return 1;
     189      if (!Chi.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiSuffix)) return 1;
     190      if (!ChiPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ChiPASSuffix)) return 1;
    173191    }
    174192  }
    … …  
    198216    prefix << dir << ShieldingPASFragmentSuffix;
    199217    if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     218    prefix.str(" ");
     219    prefix << dir << ChiFragmentSuffix;
     220    if (!ChiFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     221    prefix.str(" ");
     222    prefix << dir << ChiPASFragmentSuffix;
     223    if (!ChiPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
    200224  }
    201225
    … …  
    209233    if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
    210234    if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
     235    if (!Chi.WriteLastMatrix(argv[1], dir, ChiFragmentSuffix)) return 1;
     236    if (!ChiPAS.WriteLastMatrix(argv[1], dir, ChiPASFragmentSuffix)) return 1;
    211237  }
    212238

  • molecuilder/src/parser.hpp

    r5fd79e r348b6a  
    3131#define ShieldingFragmentSuffix ".sigma_all_fragment.all"
    3232#define ShieldingPASFragmentSuffix ".sigma_all_PAS_fragment.all"
     33#define ChiSuffix ".chi_all.csv"
     34#define ChiPASSuffix ".chi_all_PAS.csv"
     35#define ChiFragmentSuffix ".chi_all_fragment.all"
     36#define ChiPASFragmentSuffix ".chi_all_PAS_fragment.all"
    3337#define TimeSuffix ".speed"
    3438
Note: See TracChangeset for help on using the changeset viewer.