Changeset 362b0e for src/config.cpp


Ignore:
Timestamp:
Aug 18, 2008, 8:57:58 AM (17 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:
19892d
Parents:
f05407
Message:

Adaptivity fixes, MD by VerletForceIntegration introduced, MD molecule::Trajectories, atom Max::Order, no more recursive going down the fragmentation level

MD
==
molecule::Trajectories is now a map to a struct trajectory list of all the MD steps.
struct Trajectory contains STL vectors of coordinates, velocities and forces. Both are needed for the new VerletForceIntegration.
Parsing of coordinates, velocities and forces from the config file was completely rewritten:

  • in the FastParsing case, we just scan for IonType1_1 to find the last step, set the file pointer there and scan all remaining ones
  • in the other case, we create the atoms in the first step and put them in a hash for lookup on later steps and read in sequentially (with file pointer moving on).
  • This is a lot faster than the old variant.

VerletForceIntegration() implemented in a working manner with force smoothing (mean of actual and last one).
OutputTrajectoriesXYZ() now concatenates the single MD steps into one xyz file, so that the animation can be viewed with e.g. jmol or vmd
config:Deltat is now public (lazy me) and set to 1 instead of 0 initially, also Deltat is parsed accordingly (if not present, defaults to 1)
MatrixContainer::ParseMatrix from parser.cpp is used during VerletForceIntegration() to parse the forces file. Consequently, we have included parser.* in the Makefile.am.
Fix: MoleculeListClass::OutputConfigForListOfFragments() stores config file under config::configpath, before it backup'd the path twice to PathBackup.

Adaptivity
==========
Adaptivity (CheckOrderAtSite()) had compared not absolute, but real values which caused sign problems and faulty behaviour.
Adapatvity (CheckOrderAtSite()) had included atoms by Order (desired one) not by FragOrder (current one) into the list of candidates, which caused faulty behaviour.
CheckOrderAtSite() did single stepping wrong as the mask bit (last in AtomMask) was checked for true not false! Also bit was not set to false initially in FragmentMolecule().
Adaptivity: FragmentMolecule now returns 1 - continue, 2 - don't ... to tell whether we still have to continue with the adaptive cycle (is given as return value from molecuilder)
introduced atom::MaxOrder
StoreOrderAtSiteFile() now also stores the MaxOrder and ParseOrderAtSiteFromFile() parses it back into atom class

Removed Fragmentation Recursion
===============================
As we switched in the prelude of the adaptivity to a monotonous increase from order 1 to the desired one, we don't need to recursively go down each level from a given current bond order, as all these fragments have already been created on the lower orders. Consequently, FragmentBOSSANOVA does not contain NumLevels or alike anymore. This simplified the code a bit, but probably is not yet completely done. (but is working, checked on heptan).
SPFragmentGenerator() does not need Labels anymore, global ones are used for checks. Consequently, PowerSetGenerator() and FragmentSearch structure don't initialise/contain them anymore. We always compare against ...->GetTrueFather()->nr.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/config.cpp

    rf05407 r362b0e  
    4242 
    4343  MaxOuterStep=0;
    44   Deltat=0;
     44  Deltat=1;
    4545  OutVisStep=10;
    4646  OutSrcStep=5;
    … …  
    461461  int Z, No[MAX_ELEMENTS];
    462462  int verbose = 0;
    463   int MinSteps = -1;
     463  double value[3];
    464464 
    465465  /* Namen einlesen */
    … …  
    511511 
    512512  ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
    513   ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
     513  if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
     514    config::Deltat = 1;
    514515  ParseForParameter(verbose,file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    515516  ParseForParameter(verbose,file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    … …  
    643644    cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
    644645  }
    645   for (int i=0; i < config::MaxTypes; i++) {
    646     sprintf(name,"Ion_Type%i",i+1);
    647     for(int j=0;j<No[i];j++) {
    648       int repetition = 1; // which repeated keyword shall be read
    649       sprintf(keyword,"%s_%i",name, j+1);
    650       atom *neues = new atom();
    651       // then parse for each atom the coordinates as often as present
    652       if (!FastParsing) {
    653         while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
    654           ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
    655           ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional)) {
     646  int repetition = 0; // which repeated keyword shall be read
     647
     648  map<int, atom *> AtomList[config::MaxTypes];
     649  if (!FastParsing) {
     650    // parse in trajectories
     651    bool status = true;
     652    atom *neues = NULL;
     653    while (status) {
     654      for (int i=0; i < config::MaxTypes; i++) {
     655        sprintf(name,"Ion_Type%i",i+1);
     656        for(int j=0;j<No[i];j++) {
     657          sprintf(keyword,"%s_%i",name, j+1);
     658          if (repetition == 0) {
     659            neues = new atom();
     660            AtomList[i][j] = neues;
     661            neues->type = elementhash[i]; // find element type
     662            mol->AddAtom(neues);
     663          } else
     664            neues = AtomList[i][j];
     665          status = (status &&
     666                  ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
     667                  ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
     668                  ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
     669                  ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
     670          if (!status) break;
     671 
    656672          // check size of vectors
    657           if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition-1)) {
    658             cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
     673          if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
     674            //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
    659675            mol->Trajectories[neues].R.resize(repetition+10);
    660676            mol->Trajectories[neues].U.resize(repetition+10);
    661677            mol->Trajectories[neues].F.resize(repetition+10);
    662678          }
    663            
     679       
    664680          // put into trajectories list
    665681          for (int d=0;d<NDIM;d++)
    666             mol->Trajectories[neues].R.at(repetition-1).x[d] = neues->x.x[d];
     682            mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
    667683         
    668684          // parse velocities if present
    669           if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
     685          if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
    670686            neues->v.x[0] = 0.;
    671           if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
     687          if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
    672688            neues->v.x[1] = 0.;
    673           if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
     689          if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
    674690            neues->v.x[2] = 0.;
    675691          for (int d=0;d<NDIM;d++)
    676             mol->Trajectories[neues].U[repetition-1].x[d] = neues->v.x[d];
    677          
     692            mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
     693   
    678694          // parse forces if present
    679           if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &neues->v.x[0], repetition,optional))
    680             neues->v.x[0] = 0.;
    681           if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &neues->v.x[1], repetition,optional))
    682             neues->v.x[1] = 0.;
    683           if(!ParseForParameter(verbose,file, keyword, 0, 10, 1, double_type, &neues->v.x[2], repetition,optional))
    684             neues->v.x[2] = 0.;
     695          if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
     696            value[0] = 0.;
     697          if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
     698            value[1] = 0.;
     699          if(!ParseForParameter(verbose,file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
     700            value[2] = 0.;
    685701          for (int d=0;d<NDIM;d++)
    686             mol->Trajectories[neues].F[repetition-1].x[d] = neues->v.x[d];
    687          
    688           cout << "Parsed position of step " << (repetition-1) << " ";
    689           for (int d=0;d<NDIM;d++)
    690             cout << mol->Trajectories[neues].R.at(repetition-1).x[d] << "  ";          // next step
    691           cout << endl;
    692          
    693           repetition++;
     702            mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
     703
     704//            cout << "Parsed position of step " << (repetition) << ": (";
     705//            for (int d=0;d<NDIM;d++)
     706//              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
     707//            cout << ")\t(";
     708//            for (int d=0;d<NDIM;d++)
     709//              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
     710//            cout << ")\t(";
     711//            for (int d=0;d<NDIM;d++)
     712//              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
     713//            cout << ")" << endl;
    694714        }
    695         repetition--;
    696         cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
    697715      }
    698       if ((repetition < MinSteps) || (MinSteps == -1))
    699         MinSteps = repetition;
    700       ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
    701       ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
    702       ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
    703       ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
    704       if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
    705         neues->v.x[0] = 0.;
    706       if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
    707         neues->v.x[1] = 0.;
    708       if(!ParseForParameter(verbose,file, keyword, (int)FastParsing, 7, 1, double_type, &neues->v.x[2], repetition,optional))
    709         neues->v.x[2] = 0.;
    710       // here we don't care if forces are present (last in trajectories is always equal to current position)
    711       neues->type = elementhash[i]; // find element type
    712       mol->AddAtom(neues);
     716      repetition++;
    713717    }
    714   }
    715   mol->MDSteps = MinSteps;    // set equal to smallest number of parsed steps
     718    repetition--;
     719    cout << "Found " << repetition << " trajectory steps." << endl;
     720    mol->MDSteps = repetition;
     721  } else {
     722    // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
     723    repetition = 0;
     724    while ( ParseForParameter(verbose,file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
     725            ParseForParameter(verbose,file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
     726            ParseForParameter(verbose,file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
     727      repetition++;
     728    cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
     729    // parse in molecule coordinates
     730    for (int i=0; i < config::MaxTypes; i++) {
     731      sprintf(name,"Ion_Type%i",i+1);
     732      for(int j=0;j<No[i];j++) {
     733        sprintf(keyword,"%s_%i",name, j+1);
     734        atom *neues = new atom();
     735        // then parse for each atom the coordinates as often as present
     736        ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
     737        ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
     738        ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
     739        ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
     740        if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
     741          neues->v.x[0] = 0.;
     742        if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
     743          neues->v.x[1] = 0.;
     744        if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
     745          neues->v.x[2] = 0.;
     746        // here we don't care if forces are present (last in trajectories is always equal to current position)
     747        neues->type = elementhash[i]; // find element type
     748        mol->AddAtom(neues);
     749      }
     750    }
     751  }
    716752  file->close();
    717753  delete(file);
    … …  
    10151051    *output << endl;
    10161052    result = result && mol->Checkout(output);
    1017     result = result && mol->OutputTrajectories(output);
     1053    if (mol->MDSteps <=1 )
     1054      result = result && mol->Output(output);
     1055    else
     1056      result = result && mol->OutputTrajectories(output);
    10181057    return result;
    10191058  } else
    … …  
    11081147    // ... and check if it is the keyword!
    11091148    //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
    1110     if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name))) {
     1149    if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
    11111150      found++; // found the parameter!
    11121151      //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);   
    … …  
    11631202              dummy1 = dummy;
    11641203              dummy = strchr(dummy1, '\t'); // seek for tab or space
    1165               if (dummy == NULL) {
     1204              if (dummy == NULL)
    11661205                dummy = strchr(dummy1, ' ');  // if not found seek for space
    1167                 while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
    1168                   dummy++;
    1169               }
     1206              while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
     1207                dummy++;
    11701208/*              while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
    11711209                dummy++;*/
    … …  
    11891227              } else {
    11901228                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
     1229              }
     1230              if (*dummy1 == '#') {
     1231                // found comment, skipping rest of line
     1232                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
     1233                if (!sequential) { // here we need it!
     1234                  file->seekg(file_position, ios::beg);  // rewind to start position
     1235                }
     1236                Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
     1237                return 0;
    11911238              }
    11921239              //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
Note: See TracChangeset for help on using the changeset viewer.