Changeset eeec8f for src/parser.cpp


Ignore:
Timestamp:
Oct 19, 2008, 1:57:23 PM (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, 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:
84fc91
Parents:
f731ae
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/parser.cpp

    rf731ae reeec8f  
    5656MatrixContainer::MatrixContainer() {
    5757  Indices = NULL;
    58   Header = (char *) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer: *Header");
     58  Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header");
    5959  Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule
    6060  RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter");
    6161  ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter");
     62  Header[0] = NULL;
    6263  Matrix[0] = NULL;
    6364  RowCounter[0] = -1;
    … …  
    9091  Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices");
    9192 
    92   Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header");
     93  if (Header != NULL)
     94    for(int i=MatrixCounter+1;i--;)
     95      Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]");
     96  Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header");
    9397  Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    9498  Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter");
    … …  
    120124        return false;
    121125      Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]");
    122       for(int j=Matrix->RowCounter[i];j--;)
     126      for(int j=Matrix->RowCounter[i];j--;) {
    123127        Indices[i][j] = Matrix->Indices[i][j];
     128        //cout << Indices[i][j] << "\t";
     129      }
     130      //cout << endl;
    124131    }
    125132  }
    … …  
    153160    return false;
    154161  }
    155  
    156   // skip some initial lines
     162
     163  // parse header
     164  Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]");
    157165  for (int m=skiplines+1;m--;)
    158     input.getline(Header, 1023);
     166    input.getline(Header[MatrixNr], 1023);
    159167 
    160168  // scan header for number of columns
    161   line.str(Header);
     169  line.str(Header[MatrixNr]);
    162170  for(int k=skipcolumns;k--;)
    163     line >> Header;
     171    line >> Header[MatrixNr];
    164172  //cout << line.str() << endl;
    165173  ColumnCounter[MatrixNr]=0;
    … …  
    169177  }
    170178  //cout << line.str() << endl;
    171   cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
     179  //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
    172180  if (ColumnCounter[MatrixNr] == 0)
    173181    cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    … …  
    183191    }
    184192  }
    185   cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
     193  //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
    186194  if (RowCounter[MatrixNr] == 0)
    187195    cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl;
    … …  
    189197  // allocate matrix if it's not zero dimension in one direction
    190198  if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) {
    191     Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]");
     199    Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]");
    192200 
    193201    // parse in each entry for this matrix
    … …  
    195203    input.seekg(ios::beg);
    196204    for (int m=skiplines+1;m--;)
    197       input.getline(Header, 1023);    // skip header
    198     line.str(Header);
     205      input.getline(Header[MatrixNr], 1023);    // skip header
     206    line.str(Header[MatrixNr]);
    199207    for(int k=skipcolumns;k--;)  // skip columns in header too
    200208      line >> filename;
    201     strncpy(Header, line.str().c_str(), 1023); 
     209    strncpy(Header[MatrixNr], line.str().c_str(), 1023); 
    202210    for(int j=0;j<RowCounter[MatrixNr];j++) {
    203       Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     211      Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]");
    204212      input.getline(filename, 1023);
    205213      stringstream lines(filename);
    … …  
    212220      }
    213221      //cout << endl;
    214       Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseFragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");
     222      Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]");
    215223      for(int j=ColumnCounter[MatrixNr];j--;)
    216224        Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.;
    … …  
    266274 
    267275  cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl;
     276  Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule
    268277  Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    269278  RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    … …  
    271280  for(int i=MatrixCounter+1;i--;) {
    272281    Matrix[i] = NULL;
     282    Header[i] = NULL;
    273283    RowCounter[i] = -1;
    274284    ColumnCounter[i] = -1;
    … …  
    287297
    288298/** Allocates and resets the memory for a number \a MCounter of matrices.
    289  * \param *GivenHeader Header line
     299 * \param **GivenHeader Header line for each matrix
    290300 * \param MCounter number of matrices
    291301 * \param *RCounter number of rows for each matrix
    292  * \param *CCounter number of columns for each matrices
     302 * \param *CCounter number of columns for each matrix
    293303 * \return Allocation successful
    294304 */
    295 bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter)
    296 {
    297   Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader");
    298   strncpy(Header, GivenHeader, 1023);
    299 
     305bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter)
     306{
    300307  MatrixCounter = MCounter;
    301   Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule
    302   RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter");
    303   ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter");
     308  Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header");
     309  Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule
     310  RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter");
     311  ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter");
    304312  for(int i=MatrixCounter+1;i--;) {
     313    Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]");
     314    strncpy(Header[i], GivenHeader[i], 1023);
    305315    RowCounter[i] = RCounter[i];
    306316    ColumnCounter[i] = CCounter[i];
    307     Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::ParseFragmentMatrix: **Matrix[]"); 
     317    Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]"); 
    308318    for(int j=RowCounter[i]+1;j--;) {
    309       Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::ParseFragmentMatrix: *Matrix[][]");
     319      Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]");
    310320      for(int k=ColumnCounter[i];k--;)
    311321        Matrix[i][j][k] = 0.;
    … …  
    385395};
    386396
    387 /** Sums the energy with each factor and put into last element of \a ***Energies.
     397/** Sums the entries with each factor and put into last element of \a ***Matrix.
    388398 * Sums over "E"-terms to create the "F"-terms
    389399 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
    … …  
    462472      return false;
    463473    }
    464     output << Header << endl;
     474    output << Header[i] << endl;
    465475    for(int j=0;j<RowCounter[i];j++) {
    466476      for(int k=0;k<ColumnCounter[i];k++)
    … …  
    491501    return false;
    492502  }
    493   output << Header << endl;
     503  output << Header[MatrixCounter] << endl;
    494504  for(int j=0;j<RowCounter[MatrixCounter];j++) {
    495505    for(int k=0;k<ColumnCounter[MatrixCounter];k++)
    … …  
    765775      int j = Indices[ FragmentNr ][l];
    766776      if (j > RowCounter[MatrixCounter]) {
    767         cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl;
     777        cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    768778        return false;
    769779      }
    … …  
    772782          int k = Indices[ FragmentNr ][m];
    773783          if (k > ColumnCounter[MatrixCounter]) {
    774             cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << "!" << endl;
     784            cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl;
    775785            return false;
    776786          }
    … …  
    785795};
    786796
     797/** Constructor for class HessianMatrix.
     798 */
     799HessianMatrix::HessianMatrix() : MatrixContainer()
     800{
     801   IsSymmetric = true;
     802}
     803
     804/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
     805 * Sums over "E"-terms to create the "F"-terms
     806 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
     807 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order
     808 * \param Order bond order
     809 * \return true if summing was successful
     810 */
     811bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order)
     812{
     813  // go through each order
     814  for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) {
     815    //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     816    // then go per order through each suborder and pick together all the terms that contain this fragment
     817    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
     818      for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
     819        if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) {
     820          //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl;
     821          // if the fragment's indices are all in the current fragment
     822          for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
     823            int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k];
     824            //cout << "Current row index is " << k << "/" << m << "." << endl;
     825            if (m != -1) { // if it's not an added hydrogen
     826              for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
     827                //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
     828                if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) {
     829                  m = l;
     830                  break; 
     831                }
     832              }
     833              //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
     834              if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     835                cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     836                return false;
     837              }
     838             
     839              for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) {
     840                int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l];
     841                //cout << "Current column index is " << l << "/" << n << "." << endl;
     842                if (n != -1) { // if it's not an added hydrogen
     843                  for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
     844                    //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
     845                    if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) {
     846                      n = p;
     847                      break; 
     848                    }
     849                  }
     850                  //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
     851                  if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) {
     852                    cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl;
     853                    return false;
     854                  }
     855                  if (Order == SubOrder) { // equal order is always copy from Energies
     856                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     857                  } else {
     858                    Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l];
     859                  }
     860                }
     861              }
     862            }
     863            //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
     864             //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
     865          }
     866        } else {
     867          //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl;
     868        }
     869      }
     870    }
     871   //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl;
     872  }
     873 
     874  return true;
     875};
    787876
    788877/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
Note: See TracChangeset for help on using the changeset viewer.