Ignore:
Timestamp:
Jul 24, 2008, 2:00:19 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
747b10
Parents:
b3eb8e
Message:

BIG change: Hcorrection included and bugfix of analyzer (wrt forces)

  • Hcorrection is supplied via molecuilder which gives correction values if hydrogen atoms are to close and thus already feel each other's influence. This correction energy matrix is also parsed ans subtracted from the approximated energy matrix generated from the fragments. This impacts both on joiner, analyzer and molecuilder
  • Forces were not working which had multiple causes:
    • first, we stepped back down to absolute errors which rather gives a feel about convergence (1e-6 is simply small and neglegible)
    • there may have been bugs where (either in Force or in ForceFragments) adding up took place, this is cleaned up
    • more routines have been abstracted from analyzer into datacreator functions
      • the order of ions had changed, re-run with current molecuilder, calculation and subsequent joining/analyzing brought correct results finally
  • plots were bad still due to wrong axes (two more functions finding max/min values), use of wrong columns in plot file
File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/joiner.cpp

    rb3eb8e rdb3ea3  
    1919  periodentafel *periode = NULL; // and a period table of all elements
    2020  EnergyMatrix Energy;
     21  EnergyMatrix Hcorrection;
    2122  ForceMatrix Force;
    2223  EnergyMatrix EnergyFragments;
     24  EnergyMatrix HcorrectionFragments;
    2325  ForceMatrix ForceFragments;
    2426  ForceMatrix Shielding;
    … …  
    2830  KeySetsContainer KeySet; 
    2931  stringstream prefix;
     32  char *dir = NULL;
     33  bool Hcorrected = true;
    3034
    3135  cout << "Joiner" << endl;
    … …  
    3943    cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl;
    4044    return 1;
     45  } else {
     46    dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir");
     47    strcpy(dir, "/");
     48    strcat(dir, argv[2]);
    4149  }
    4250  if (argc > 3) {
    … …  
    5260 
    5361  // ------------- Parse through all Fragment subdirs --------
    54   if (!Energy.ParseMatrix(argv[1], argv[2], EnergySuffix, 0,0)) return 1;
    55   if (!Force.ParseMatrix(argv[1], argv[2], ForcesSuffix, 0,0)) return 1;
     62  if (!Energy.ParseMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1;
     63  Hcorrected = Hcorrection.ParseMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0);
     64  if (!Force.ParseMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1;
    5665  if (periode != NULL) { // also look for PAS values
    57     if (!Shielding.ParseMatrix(argv[1], argv[2], ShieldingSuffix, 1, 0)) return 1;
    58     if (!ShieldingPAS.ParseMatrix(argv[1], argv[2], ShieldingPASSuffix, 1, 0)) return 1;
     66    if (!Shielding.ParseMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1;
     67    if (!ShieldingPAS.ParseMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1;
    5968  }
    6069
    6170  // ---------- Parse the TE Factors into an array -----------------
    6271  if (!Energy.ParseIndices()) return 1;
     72  if (Hcorrected) Hcorrection.ParseIndices();
    6373 
    6474  // ---------- Parse the Force indices into an array ---------------
    … …  
    7686  if (!KeySet.ParseManyBodyTerms()) return 1;
    7787  if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1;
     88  if (Hcorrected)  HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter);
    7889  if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1;
    7990  if (periode != NULL) { // also look for PAS values
    … …  
    8697  if(!Force.SetLastMatrix(0., 2)) return 1;
    8798  if (periode != NULL) { // also look for PAS values
    88     if(!Shielding.SetLastMatrix(0., 0)) return 1;
     99    if(!Shielding.SetLastMatrix(0., 2)) return 1;
    89100    if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1;
    90101  }
    … …  
    97108    cout << "Summing energy of order " << BondOrder+1 << " ..." << endl;
    98109    if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1;
    99     if (!Energy.SumSubEnergy(EnergyFragments, KeySet, BondOrder)) return 1;
     110    if (Hcorrected) {
     111      HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder);
     112      if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1;
     113      if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);
     114    } else
     115      if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1;
    100116    // --------- sum up Forces --------------------
    101117    cout << "Summing forces of order " << BondOrder+1 << " ..." << endl;
    102118    if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1;
    103     if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder)) return 1;
     119    if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1;
    104120    if (periode != NULL) { // also look for PAS values
    105121      cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl;
    106122      if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1;
    107       if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder)) return 1;
     123      if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1;
    108124      if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1;
    109       if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder)) return 1;
     125      if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1;
    110126    }
    111127
    112128    // --------- write the energy and forces file --------------------
    113129    prefix.str(" ");
    114     prefix << argv[2] << OrderSuffix << (BondOrder+1);
     130    prefix << dir << OrderSuffix << (BondOrder+1);
    115131    cout << "Writing files " << argv[1] << prefix.str() << ". ..." << endl;
    116132    // energy
    … …  
    126142  // fragments
    127143  prefix.str(" ");
    128   prefix << argv[2] << EnergyFragmentSuffix;
     144  prefix << dir << EnergyFragmentSuffix;
    129145  if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     146  if (Hcorrected) {
     147    prefix.str(" ");
     148    prefix << dir << HcorrectionFragmentSuffix;
     149    if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
     150  }
    130151  prefix.str(" ");
    131   prefix << argv[2] << ForceFragmentSuffix;
     152  prefix << dir << ForceFragmentSuffix;
    132153  if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
    133154  if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], FRAGMENTPREFIX ENERGYPERFRAGMENT, "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1;
    134155  if (periode != NULL) { // also look for PAS values
    135156    prefix.str(" ");
    136     prefix << argv[2] << ShieldingFragmentSuffix;
     157    prefix << dir << ShieldingFragmentSuffix;
    137158    if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
    138159    prefix.str(" ");
    139     prefix << argv[2] << ShieldingPASFragmentSuffix;
     160    prefix << dir << ShieldingPASFragmentSuffix;
    140161    if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1;
    141162  }
    142163
    143164  // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds
    144   if (!Energy.WriteLastMatrix(argv[1], argv[2], EnergyFragmentSuffix)) return 1;
    145   if (!Force.WriteLastMatrix(argv[1], argv[2], ForceFragmentSuffix)) return 1;
     165  if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1;
     166  if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);
     167  if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1;
    146168  if (periode != NULL) { // also look for PAS values
    147     if (!Shielding.WriteLastMatrix(argv[1], argv[2], ShieldingFragmentSuffix)) return 1;
    148     if (!ShieldingPAS.WriteLastMatrix(argv[1], argv[2], ShieldingPASFragmentSuffix)) return 1;
     169    if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1;
     170    if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1;
    149171  }
    150172
    151173  // exit 
    152174  delete(periode);
     175  Free((void **)&dir, "main: *dir");
    153176  cout << "done." << endl;
    154177  return 0;
Note: See TracChangeset for help on using the changeset viewer.