Changeset 362b0e


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, 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.

Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rf05407 r362b0e  
    1 SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
    2 HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp periodentafel.hpp stackclass.hpp vector.hpp
     1SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp parser.cpp periodentafel.cpp vector.cpp verbose.cpp
     2HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp vector.hpp
    33
    44bin_PROGRAMS = molecuilder joiner analyzer

  • src/atom.cpp

    rf05407 r362b0e  
    2828  LowpointNr = -1;
    2929  AdaptiveOrder = 0;
     30  MaxOrder = false;
    3031};
    3132

  • src/builder.cpp

    rf05407 r362b0e  
    703703    output.open(filename, ios::trunc);
    704704  }
    705   if (mol->OutputXYZ(&output))
    706     cout << Verbose(0) << "Saving of XYZ file successful." << endl;
    707   else
    708     cout << Verbose(0) << "Saving of XYZ file failed." << endl;
     705  if (mol->MDSteps <= 1) {
     706    if (mol->OutputXYZ(&output))
     707      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
     708    else
     709      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
     710  } else {
     711    if (mol->OutputTrajectoriesXYZ(&output))
     712      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
     713    else
     714      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
     715  }
    709716  output.close();
    710717  output.clear();
    … …  
    733740  string line;
    734741  atom *first;
    735   int ExitFlag = 0;
     742  bool SaveFlag = false;
     743  int ExitFlag = 1;
    736744  int j;
    737745  double volume = 0.;
    … …  
    758766            cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
    759767            cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
     768            cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    760769            cout << "\t-O\tCenter atoms in origin." << endl;
    761770            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    762771            cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    763             cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl;
     772            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
    764773            cout << "\t-h/-H/-?\tGive this help screen." << endl;
    765774            cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    767776            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    768777            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
    769             cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl;
     778            cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
    770779            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    771780            cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
    … …  
    860869          switch(argv[argptr-1][1]) {
    861870            case 'p':
    862               ExitFlag = 1;
     871              SaveFlag = true;
    863872              cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
    864873              if (!mol->AddXYZFile(argv[argptr]))
    … …  
    875884        if (config_present == present) {
    876885          switch(argv[argptr-1][1]) {
     886            case 'D':
     887              {
     888                cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
     889                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
     890                int *MinimumRingSize = NULL;
     891                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
     892                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
     893                Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
     894                delete[](MinimumRingSize);
     895                if (Subgraphs != NULL) {
     896                  while (Subgraphs->next != NULL) {
     897                    Subgraphs = Subgraphs->next;
     898                    delete(Subgraphs->previous);
     899                  }
     900                  delete(Subgraphs);
     901                }
     902              }
     903              argptr+=1;
     904              break;
    877905            case 'P':
    878               ExitFlag = 1;
     906              SaveFlag = true;
    879907              cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
    880               if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1])))
     908              if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
    881909                cout << Verbose(2) << "File not found." << endl;
    882910              else
    883911                cout << Verbose(2) << "File found and parsed." << endl;
    884               argptr+=2;
     912              argptr+=1;
    885913              break;
    886914            case 't':
    887               ExitFlag = 1;
     915              SaveFlag = true;
    888916              cout << Verbose(1) << "Translating all ions to new origin." << endl;
    889917              for (int i=NDIM;i--;)
    … …  
    893921              break;
    894922            case 'a':
    895               ExitFlag = 1;
     923              SaveFlag = true;
    896924              cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
    897925              first = new atom;
    … …  
    908936              break;
    909937            case 's':
    910               ExitFlag = 1;
     938              SaveFlag = true;
    911939              j = -1;
    912940              cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
    … …  
    929957              break;
    930958            case 'b':
    931               ExitFlag = 1;
     959              SaveFlag = true;
    932960              j = -1;
    933961              cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
    … …  
    944972              break;
    945973            case 'c':
    946               ExitFlag = 1;
     974              SaveFlag = true;
    947975              j = -1;
    948976              cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
    … …  
    961989              break;
    962990            case 'O':
    963               ExitFlag = 1;
     991              SaveFlag = true;
    964992              cout << Verbose(1) << "Centering atoms in origin." << endl;
    965993              mol->CenterOrigin((ofstream *)&cout, &x);
    … …  
    967995              break;
    968996            case 'r':
    969               ExitFlag = 1;
     997              SaveFlag = true;
    970998              cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
    971999              break;
    9721000            case 'F':
    9731001            case 'f':
    974               if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
    9751002              cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
    9761003              if (argc >= argptr+2) {
    … …  
    9801007                cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
    9811008                if (mol->first->next != mol->last) {
    982                   mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
     1009                  ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
    9831010                }
    9841011                end = clock();
    … …  
    9901017              break;
    9911018            case 'm':
    992               ExitFlag = 1;
    9931019              j = atoi(argv[argptr++]);
    9941020              if ((j<0) || (j>1)) {
    … …  
    9961022                j = 0;
    9971023              }
    998               if (j)
     1024              if (j) {
     1025                SaveFlag = true;
    9991026                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
    1000               else
     1027              } else
    10011028                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
    10021029              mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
    10031030              break;
    10041031            case 'o':
    1005               ExitFlag = 1;
     1032              SaveFlag = true;
    10061033              cout << Verbose(0) << "Evaluating volume of the convex envelope.";
    10071034              VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
    … …  
    10131040              {
    10141041                double density;
    1015                 ExitFlag = 1;
     1042                SaveFlag = true;
    10161043                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
    10171044                density = atof(argv[argptr++]);
    … …  
    10301057              break;
    10311058            case 'd':
     1059              SaveFlag = true;
    10321060              for (int axis = 1; axis <= NDIM; axis++) {
    10331061                int faktor = atoi(argv[argptr++]);
    … …  
    10891117      } else argptr++;
    10901118    } while (argptr < argc);
    1091     if (ExitFlag == 1) // 1 means save and exit
     1119    if (SaveFlag)
    10921120      SaveConfig(ConfigFileName, &configuration, periode, mol);
    1093     if (ExitFlag >= 1) { // 2 means just exit
     1121    if ((ExitFlag >= 1)) {
    10941122      delete(mol);
    10951123      delete(periode);
    1096       return (1);
     1124      return (ExitFlag);
    10971125    }
    10981126  } else {  // no arguments, hence scan the elements db

  • 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);

  • src/moleculelist.cpp

    rf05407 r362b0e  
    427427    *out << endl;
    428428    // prepare output of config file
    429     sprintf(FragmentName, "%s/%s%s.conf", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     429    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    430430    outputFragment.open(FragmentName, ios::out);
    431431    //strcpy(PathBackup, configuration->configpath);
    … …  
    739739void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
    740740{
     741  *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
    741742  KeySet *TempSet = new KeySet;
    742743  if (FragmentList[FragmentCounter] != NULL) {
    … …  
    754755    next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
    755756  FragmentCounter--;
     757  *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
    756758};
    757759

  • src/molecules.cpp

    rf05407 r362b0e  
    455455  getline(xyzfile,line,'\n'); // Read comment
    456456  cout << Verbose(1) << "Comment: " << line << endl;
    457 
     457 
     458  if (MDSteps == 0) // no atoms yet present
     459    MDSteps++;
    458460  for(i=0;i<NumberOfAtoms;i++){
    459461    Walker = new atom;
    … …  
    471473      Walker->type = elemente->FindElement(1);
    472474    }
    473     for(j=NDIM;j--;)
     475    if (Trajectories[Walker].R.size() <= (unsigned int)MDSteps) {
     476      Trajectories[Walker].R.resize(MDSteps+10);
     477      Trajectories[Walker].U.resize(MDSteps+10);
     478      Trajectories[Walker].F.resize(MDSteps+10);
     479    }
     480    for(j=NDIM;j--;) {
    474481      Walker->x.x[j] = x[j];
    475      AddAtom(Walker);  // add to molecule
     482      Trajectories[Walker].R.at(MDSteps-1).x[j] = x[j];
     483      Trajectories[Walker].U.at(MDSteps-1).x[j] = 0;
     484      Trajectories[Walker].F.at(MDSteps-1).x[j] = 0;
     485    }
     486    AddAtom(Walker);  // add to molecule
    476487    delete(item);
    477488  }
    … …  
    965976
    966977/** Parses nuclear forces from file and performs Verlet integration.
     978 * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
     979 * have to transform them).
    967980 * This adds a new MD step to the config file.
    968981 * \param *file filename
    969982 * \param delta_t time step width in atomic units
     983 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    970984 * \return true - file found and parsed, false - file not found or imparsable
    971985 */
    972 bool molecule::VerletForceIntegration(char *file, double delta_t)
    973 {
    974   bool status = true;
     986bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
     987{
    975988  element *runner = elemente->start;
    976989  atom *walker = NULL;
    977   int ElementNo, AtomNo;
     990  int AtomNo;
    978991  ifstream input(file);
    979   double Forcesvector[3];
    980   double dummy;
    981992  string token;
    982   size_t oldpos, newpos;
    983993  stringstream item;
    984   int i, Ion[2];
    985   double a, IonMass;
    986   unsigned int size;
     994  double a, IonMass, Vector[NDIM];
     995  ForceMatrix Force;
    987996
    988997  CountElements();  // make sure ElementsInMolecule is up to date
    … …  
    9921001    return false;
    9931002  } else {
    994     ElementNo = 0;
     1003    // parse file into ForceMatrix
     1004    if (!Force.ParseMatrix(file, 0,0,0)) {
     1005      cerr << "Could not parse Force Matrix file " << file << "." << endl;
     1006      return false;
     1007    }
     1008    if (Force.RowCounter[0] != AtomCount) {
     1009      cerr << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl;
     1010      return false;
     1011    }
     1012    // correct Forces
     1013    for(int d=0;d<NDIM;d++)
     1014      Vector[d] = 0.;
     1015    for(int i=0;i<AtomCount;i++)
     1016      for(int d=0;d<NDIM;d++) {
     1017        Vector[d] += Force.Matrix[0][i][d+5];
     1018      }
     1019    for(int i=0;i<AtomCount;i++)
     1020      for(int d=0;d<NDIM;d++) {
     1021        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
     1022      }
     1023    // and perform Verlet integration for each atom with position, velocity and force vector
     1024    runner = elemente->start;
    9951025    while (runner->next != elemente->end) { // go through every element
    9961026      runner = runner->next;
    9971027      IonMass = runner->mass;
     1028      a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    9981029      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
    999         ElementNo++;
    10001030        AtomNo = 0;
    10011031        walker = start;
    … …  
    10031033          walker = walker->next;
    10041034          if (walker->type == runner) { // if this atom fits to element
    1005             AtomNo++;
    1006             // parse Forcevector for this ion
    1007             getline (input, token, '\n');
    1008             newpos = oldpos = i = 0;
    1009             while ((newpos = token.find( '\t', oldpos)) != oldpos) {
    1010               i++;
    1011               item.str();
    1012               switch (i) {
    1013                 case 1:
    1014                 case 2:
    1015                   item >> Ion[i-1];
    1016                   break;
    1017                 case 6:
    1018                 case 7:
    1019                 case 8:
    1020                   item >> Forcesvector[i-6];
    1021                   break;
    1022                 case 3:
    1023                   if ((Ion[0] != ElementNo) || (Ion[1] != AtomNo))
    1024                     cout << "ERROR: Expected " << ElementNo << ", " << AtomNo << "but parsed " << Ion[0] << ", " << Ion[1] << "." << endl;
    1025                   // still parse into dummy! Hence no break here ...
    1026                 default:
    1027                   item >> dummy;
    1028               }
    1029               oldpos = newpos;
     1035            // check size of vectors
     1036            if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
     1037              //cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
     1038              Trajectories[walker].R.resize(MDSteps+10);
     1039              Trajectories[walker].U.resize(MDSteps+10);
     1040              Trajectories[walker].F.resize(MDSteps+10);
    10301041            }
    10311042           
    1032             // perform Verlet integration for this atom with position, velocity and force vector
    1033            
    1034             // UpdateR
    1035 //            for (int d=0; d<NDIM; d++) {
    1036 //              R_old_old[d] = R_old[d];        // shift old values
    1037 //              R_old[d] = R[d];
    1038 //              R[d] += delta_t*(U[d]+a*FIon[d]);     // F = m * a and s = 0.5 * a * t^2
    1039 //              FIon_old[d] = FIon[d];          // store old force
    1040 //            }
    1041 //            // UpdateU
    1042 //            a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    1043 //            for (int d=0; d<NDIM; d++) {
    1044 //                U[d] += a * (FIon[d]+FIon_old[d]);
    1045 //            }
    1046            
    1047            
    1048             // check size of vectors
    1049             size = Trajectories[walker].R.size();
    1050             if (size <= (unsigned int)(MDSteps)) {
    1051               cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    1052               Trajectories[walker].R.resize(size+10);
    1053               Trajectories[walker].U.resize(size+10);
    1054               Trajectories[walker].F.resize(size+10);
     1043            // Update R (and F)
     1044            for (int d=0; d<NDIM; d++) {
     1045              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
     1046              Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
     1047              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
     1048              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
    10551049            }
    1056             // add trajectory point
    1057 
     1050            // Update U
     1051            for (int d=0; d<NDIM; d++) {
     1052              Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
     1053              Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
     1054            }
     1055//            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
     1056//            for (int d=0;d<NDIM;d++)
     1057//              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
     1058//            cout << ")\t(";
     1059//            for (int d=0;d<NDIM;d++)
     1060//              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
     1061//            cout << ")" << endl;
     1062            // next atom
     1063            AtomNo++;
    10581064          }
    10591065        }
    10601066      }
    10611067    }
    1062     return true;
    1063   }
     1068  }
     1069//  // correct velocities (rather momenta) so that center of mass remains motionless
     1070//  for(int d=0;d<NDIM;d++)
     1071//    Vector[d] = 0.;
     1072//  IonMass = 0.;
     1073//  walker = start;
     1074//  while (walker->next != end) { // go through every atom
     1075//    walker = walker->next;
     1076//    IonMass += walker->type->mass;  // sum up total mass
     1077//    for(int d=0;d<NDIM;d++) {
     1078//      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
     1079//    }
     1080//  }
     1081//  walker = start;
     1082//  while (walker->next != end) { // go through every atom of this element
     1083//    walker = walker->next;
     1084//    for(int d=0;d<NDIM;d++) {
     1085//      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
     1086//    }
     1087//  }
     1088  MDSteps++;
     1089
    10641090 
    10651091  // exit
    1066   return status;
     1092  return true;
    10671093};
    10681094
    … …  
    13721398              if (Trajectories[walker].F.at(step).Norm() > MYEPSILON)
    13731399                *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t";
    1374               *out << " # Number in molecule " << walker->nr << endl;
     1400              *out << "\t# Number in molecule " << walker->nr << endl;
    13751401            }
    13761402          }
    … …  
    14131439};
    14141440
     1441/** Prints molecule with all its trajectories to *out as xyz file.
     1442 * \param *out output stream
     1443 */
     1444bool molecule::OutputTrajectoriesXYZ(ofstream *out)
     1445{
     1446  atom *walker = NULL;
     1447  int No = 0;
     1448  time_t now;
     1449 
     1450  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     1451  walker = start;
     1452  while (walker->next != end) { // go through every atom and count
     1453    walker = walker->next;
     1454    No++;
     1455  }
     1456  if (out != NULL) {
     1457    for (int step=0;step<MDSteps;step++) {
     1458      *out << No << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
     1459      walker = start;
     1460      while (walker->next != end) { // go through every atom of this element
     1461        walker = walker->next;
     1462        *out << walker->type->symbol << "\t" << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2] << endl;
     1463      }
     1464    }
     1465    return true;
     1466  } else
     1467    return false;
     1468};
     1469
    14151470/** Prints molecule to *out as xyz file.
    1416  * \param *out output stream
     1471* \param *out output stream
    14171472 */
    14181473bool molecule::OutputXYZ(ofstream *out) const
    … …  
    26152670        lines++;
    26162671      }
    2617       *out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
     2672      //*out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
    26182673      InputFile.clear();
    26192674      InputFile.seekg(ios::beg);
    … …  
    26272682        InputFile.getline(buffer, MAXSTRINGSIZE);
    26282683        if (strlen(buffer) > 2) {
    2629           //*out << Verbose(2) << "Scanning: " << buffer;
     2684          //*out << Verbose(2) << "Scanning: " << buffer << endl;
    26302685          stringstream line(buffer);
    26312686          line >> FragOrder;
    … …  
    26342689          line >> ws >> Value;
    26352690          No -= 1;  // indices start at 1 in file, not 0
    2636           //*out << Verbose(2) << " - yields (" << No << "," << Value << ")" << endl;
     2691          //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
    26372692
    26382693          // clean the list of those entries that have been superceded by higher order terms already
    26392694          map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
    26402695          if (marker != IndexKeySetList.end()) {  // if found
     2696            Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes em not equal without changing anything actually
    26412697            // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
    2642             pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( Value, Order) ));
     2698            pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
    26432699            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
    26442700            if (!InsertedElement.second) { // this root is already present
    … …  
    26462702                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
    26472703                {  // if value is smaller, update value and order
    2648                 (*PresentItem).second.first = Value;
     2704                (*PresentItem).second.first = fabs(Value);
    26492705                (*PresentItem).second.second = FragOrder;
    26502706                *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
     2707              } else {
     2708                *out << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl;
    26512709              }
    26522710            } else {
    … …  
    26642722        Walker = FindAtom((*runner).first);
    26652723        if (Walker != NULL) {
    2666           if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
     2724          //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
     2725          if (!Walker->MaxOrder) {
    26672726            *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
    26682727            FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
     2728          } else {
     2729            *out << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl;
    26692730          }
    26702731        } else {
    … …  
    27132774      }
    27142775    }
    2715     if ((Order == 0) && (AtomMask[AtomCount] == true))  // single stepping, just check
     2776    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
    27162777      status = true;
    27172778     
    2718     if (!status)
    2719       *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
     2779    if (!status) {
     2780      if (Order == 0)
     2781        *out << Verbose(1) << "Single stepping done." << endl;
     2782      else
     2783        *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
     2784    }
    27202785  }
    27212786 
    … …  
    27322797};
    27332798
    2734 /** Create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file.
     2799/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
    27352800 * \param *out output stream for debugging
    27362801 * \param *&SortIndex Mapping array of size molecule::AtomCount
    … …  
    27812846 * \param *configuration configuration for writing config files for each fragment
    27822847 */
    2783 void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
     2848int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
    27842849{
    27852850        MoleculeListClass *BondFragments = NULL;
    … …  
    27912856  fstream File;
    27922857  bool FragmentationToDo = true;
     2858  bool CheckOrder = false;
    27932859  Graph **FragmentList = NULL;
    27942860  Graph *ParsedFragmentList = NULL;
    … …  
    28372903  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    28382904  AtomMask = new bool[AtomCount+1];
    2839   while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath)) {
     2905  AtomMask[AtomCount] = false;
     2906  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
     2907  while ((CheckOrder = CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
     2908    FragmentationToDo = FragmentationToDo || CheckOrder;
    28402909    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    28412910    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    … …  
    28892958
    28902959  // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
    2891   // allocate memory for the pointer array and transmorph graphs into full molecular fragments
    2892   BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
    2893   int k=0;
    2894   for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    2895     KeySet test = (*runner).first;
    2896     *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
    2897     BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
    2898     k++;
    2899   }
    2900   *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
    2901  
    2902   // ===== 9. Save fragments' configuration and keyset files et al to disk ===
    2903   if (BondFragments->NumberOfMolecules != 0) {
    2904     // create the SortIndex from BFS labels to order in the config file
    2905     CreateMappingLabelsToConfigSequence(out, SortIndex);
     2960  //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
     2961    // allocate memory for the pointer array and transmorph graphs into full molecular fragments
     2962    BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
     2963    int k=0;
     2964    for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
     2965      KeySet test = (*runner).first;
     2966      *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
     2967      BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
     2968      k++;
     2969    }
     2970    *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
    29062971   
    2907     *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
    2908     if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
    2909       *out << Verbose(1) << "All configs written." << endl;
    2910     else
    2911       *out << Verbose(1) << "Some config writing failed." << endl;
    2912 
    2913     // store force index reference file
    2914     BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
    2915    
    2916     // store keysets file
    2917     StoreKeySetFile(out, TotalGraph, configuration->configpath);
    2918 
    2919     // store Adjacency file
    2920     StoreAdjacencyToFile(out, configuration->configpath);
    2921 
    2922     // store Hydrogen saturation correction file
    2923     BondFragments->AddHydrogenCorrection(out, configuration->configpath);
    2924    
    2925     // store adaptive orders into file
    2926     StoreOrderAtSiteFile(out, configuration->configpath);
    2927    
    2928     // restore orbital and Stop values
    2929     CalculateOrbitals(*configuration);
    2930    
    2931     // free memory for bond part
    2932     *out << Verbose(1) << "Freeing bond memory" << endl;
    2933     delete(FragmentList); // remove bond molecule from memory
    2934     Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
    2935   } else
    2936     *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
    2937  
     2972    // ===== 9. Save fragments' configuration and keyset files et al to disk ===
     2973    if (BondFragments->NumberOfMolecules != 0) {
     2974      // create the SortIndex from BFS labels to order in the config file
     2975      CreateMappingLabelsToConfigSequence(out, SortIndex);
     2976     
     2977      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
     2978      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
     2979        *out << Verbose(1) << "All configs written." << endl;
     2980      else
     2981        *out << Verbose(1) << "Some config writing failed." << endl;
     2982 
     2983      // store force index reference file
     2984      BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
     2985     
     2986      // store keysets file
     2987      StoreKeySetFile(out, TotalGraph, configuration->configpath);
     2988 
     2989      // store Adjacency file
     2990      StoreAdjacencyToFile(out, configuration->configpath);
     2991 
     2992      // store Hydrogen saturation correction file
     2993      BondFragments->AddHydrogenCorrection(out, configuration->configpath);
     2994     
     2995      // store adaptive orders into file
     2996      StoreOrderAtSiteFile(out, configuration->configpath);
     2997     
     2998      // restore orbital and Stop values
     2999      CalculateOrbitals(*configuration);
     3000     
     3001      // free memory for bond part
     3002      *out << Verbose(1) << "Freeing bond memory" << endl;
     3003      delete(FragmentList); // remove bond molecule from memory
     3004      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
     3005    } else
     3006      *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
     3007  //} else
     3008  //  *out << Verbose(1) << "No fragments to store." << endl;
    29383009  *out << Verbose(0) << "End of bond fragmentation." << endl;
     3010
     3011  return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
    29393012};
    29403013
    … …  
    29573030    while (Walker->next != end) {
    29583031      Walker = Walker->next;
    2959       file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << endl;
    2960       *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "." << endl;
     3032      file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << endl;
     3033      *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << "." << endl;
    29613034    }
    29623035    file.close();
    … …  
    29793052{
    29803053  unsigned char *OrderArray = (unsigned char *) Malloc(sizeof(unsigned char)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
     3054  bool *MaxArray = (bool *) Malloc(sizeof(bool)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
    29813055  bool status;
    2982   int AtomNr;
     3056  int AtomNr, value;
    29833057  stringstream line;
    29843058  ifstream file;
    2985   int Order;
    29863059
    29873060  *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
    … …  
    29913064  file.open(line.str().c_str());
    29923065  if (file != NULL) {
    2993     for (int i=AtomCount;i--;) // initialise with 0
     3066    for (int i=AtomCount;i--;) { // initialise with 0
    29943067      OrderArray[i] = 0;
     3068      MaxArray[i] = 0;
     3069    }
    29953070    while (!file.eof()) { // parse from file
     3071      AtomNr = -1;
    29963072      file >> AtomNr;
    2997       file >> Order;
    2998       OrderArray[AtomNr] = (unsigned char) Order;
    2999       //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << "." << endl;
     3073      if (AtomNr != -1) {   // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
     3074        file >> value;
     3075        OrderArray[AtomNr] = value;
     3076        file >> value;
     3077        MaxArray[AtomNr] = value;
     3078        //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
     3079      }
    30003080    }
    30013081    atom *Walker = start;
    … …  
    30033083      Walker = Walker->next;
    30043084      Walker->AdaptiveOrder = OrderArray[Walker->nr];
    3005       *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << "." << endl;
     3085      Walker->MaxOrder = MaxArray[Walker->nr];
     3086      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
    30063087    }
    30073088    file.close();
    … …  
    30133094  }
    30143095  Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
     3096  Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
    30153097 
    30163098  *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
    … …  
    36023684  int FragmentCounter;
    36033685  int CurrentIndex;
    3604   int *Labels;
    36053686  double TEFactor;
    36063687  int *ShortestPathList;
    … …  
    36703751                OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add
    36713752          //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
    3672           //if ((!FragmentSearch->UsedList[OtherWalker->nr][i]) && (FragmentSearch->Labels[OtherWalker->nr] > FragmentSearch->Labels[FragmentSearch->Root->nr])) {
    3673             //*out << Verbose(2+verbosity) << "Not used so far, label " << FragmentSearch->Labels[OtherWalker->nr] << " is bigger than Root's " << FragmentSearch->Labels[FragmentSearch->Root->nr] << "." << endl;
    36743753          *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
    36753754          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
    … …  
    37303809        InsertFragmentIntoGraph(out, FragmentSearch);
    37313810        //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
    3732         //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList,FragmentSearch->Labels, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
     3811        //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
    37333812      }
    37343813     
    … …  
    38213900int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
    38223901{
    3823   int SP, UniqueIndex, AtomKeyNr;
     3902  int SP, AtomKeyNr;
    38243903  atom *Walker = NULL, *OtherWalker = NULL, *Predecessor = NULL;
    38253904  bond *Binder = NULL;
    38263905  bond *CurrentEdge = NULL;
    38273906  bond **BondsList = NULL;
    3828   int RootKeyNr = FragmentSearch.Root->nr;
     3907  int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
    38293908  int Counter = FragmentSearch.FragmentCounter;
    38303909  int RemainingWalkers;
    … …  
    38343913
    38353914  // prepare Label and SP arrays of the BFS search
    3836   UniqueIndex = 0;
    3837   if (FragmentSearch.Labels[RootKeyNr] == -1)
    3838     FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
    3839   FragmentSearch.ShortestPathList[RootKeyNr] = 0;
     3915  FragmentSearch.ShortestPathList[FragmentSearch.Root->nr] = 0;
    38403916
    38413917  // prepare root level (SP = 0) and a loop bond denoting Root
    … …  
    38713947      Predecessor = CurrentEdge->leftatom;    // ... and leftatom is predecessor
    38723948      AtomKeyNr = Walker->nr;
    3873       *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
     3949      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
    38743950      // check for new sp level
    38753951      // go through all its bonds
    … …  
    38853961          *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
    38863962          // set the label if not set (and push on root stack as well)
    3887           if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
    3888             FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
    3889             *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
    3890           } else {
    3891             *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
    3892           }
    3893           if ((OtherWalker != Predecessor) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
     3963          if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
    38943964            FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
    38953965            *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
    … …  
    38993969            FragmentSearch.BondsPerSPCount[SP+1]++;
    39003970            *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl;
    3901           } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor " << *Predecessor << "." << endl;
     3971          } else {
     3972            if (OtherWalker != Predecessor)
     3973              *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl;
     3974            else
     3975              *out << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl;
     3976          }
    39023977        } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
    39033978      }
    … …  
    40034078      for (int i=NDIM;i--;) {
    40044079        tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
    4005         //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
     4080        *out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
    40064081        if (tmp > BondDistance) {
    40074082          OtherBinder = Binder->next; // note down binding partner for later re-insertion
    40084083          unlink(Binder);   // unlink bond
    4009           //*out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
     4084          *out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
    40104085          flag = true;
    40114086          break;
    … …  
    41954270  // initialise the fragments structure
    41964271  FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
    4197   FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
    41984272  FragmentSearch.FragmentCounter = 0;
    41994273  FragmentSearch.FragmentSet = new KeySet;
    42004274  FragmentSearch.Root = FindAtom(RootKeyNr);
    42014275  for (int i=AtomCount;i--;) {
    4202     FragmentSearch.Labels[i] = -1;
    42034276    FragmentSearch.ShortestPathList[i] = -1;
    42044277  }
    … …  
    42494322      // create top order where nothing is reduced
    42504323      *out << Verbose(0) << "==============================================================================================================" << endl;
    4251       *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl;
     4324      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "
    42524325 
    42534326      // Create list of Graphs of current Bond Order (i.e. F_{ij})
    … …  
    42604333      if (NumMoleculesOfOrder[RootNr] != 0) {
    42614334        NumMolecules = 0;
     4335
     4336        // we don't have to dive into suborders! These keysets are all already created on lower orders!
     4337        // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)
    42624338       
    4263         if ((NumLevels >> 1) > 0) {
    4264           // create lower order fragments
    4265           *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
    4266           Order = Walker->AdaptiveOrder;
    4267           for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
    4268             // step down to next order at (virtual) boundary of powers of 2 in array
    4269             while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
    4270               Order--;
    4271             *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
    4272             for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
    4273               int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
    4274               *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
    4275               *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
    4276        
    4277               // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
    4278               //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
    4279               //NumMolecules = 0;
    4280               FragmentLowerOrdersList[RootNr][dest] = new Graph;
    4281               for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
    4282                 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    4283                   Graph TempFragmentList;
    4284                   FragmentSearch.TEFactor = -(*runner).second.second;
    4285                   FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
    4286                   FragmentSearch.Root = FindAtom(*sprinter);
    4287                   NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
    4288                   // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
    4289                   *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
    4290                   InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
    4291                 }
    4292               }
    4293               *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
    4294             }
    4295           }
    4296         }
     4339//        if ((NumLevels >> 1) > 0) {
     4340//          // create lower order fragments
     4341//          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
     4342//          Order = Walker->AdaptiveOrder;
     4343//          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
     4344//            // step down to next order at (virtual) boundary of powers of 2 in array
     4345//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))   
     4346//              Order--;
     4347//            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
     4348//            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
     4349//              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
     4350//              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
     4351//              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
     4352//       
     4353//              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
     4354//              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
     4355//              //NumMolecules = 0;
     4356//              FragmentLowerOrdersList[RootNr][dest] = new Graph;
     4357//              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
     4358//                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
     4359//                  Graph TempFragmentList;
     4360//                  FragmentSearch.TEFactor = -(*runner).second.second;
     4361//                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
     4362//                  FragmentSearch.Root = FindAtom(*sprinter);
     4363//                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
     4364//                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
     4365//                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
     4366//                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
     4367//                }
     4368//              }
     4369//              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
     4370//            }
     4371//          }
     4372//        }
    42974373      } else {
    4298         *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
     4374        Walker->GetTrueFather()->MaxOrder = true;
     4375//        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
    42994376      }
    43004377      // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
    43014378      //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
    43024379      TotalNumMolecules += NumMoleculesOfOrder[RootNr];
    4303       *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
     4380//      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
    43044381      RootStack.push_back(RootKeyNr); // put back on stack
    43054382      RootNr++;
    … …  
    43204397  // cleanup FragmentSearch structure
    43214398  Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
    4322   Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
    43234399  delete(FragmentSearch.FragmentSet);
    43244400 

  • src/molecules.hpp

    rf05407 r362b0e  
    2424
    2525#include "helpers.hpp"
     26#include "parser.hpp"
    2627#include "periodentafel.hpp"
    2728#include "stackclass.hpp"
    … …  
    113114    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
    114115    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
     116    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
    115117 
    116118  atom();
    … …  
    231233        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
    232234        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
    233         bool VerletForceIntegration(char *file, double delta_t);
     235        bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
    234236       
    235237  bool CheckBounds(const Vector *x) const;
    … …  
    256258 
    257259  /// Fragment molecule by two different approaches:
    258   void FragmentMolecule(ofstream *out, int Order, config *configuration);
     260  int FragmentMolecule(ofstream *out, int Order, config *configuration);
    259261  bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
    260262  bool StoreAdjacencyToFile(ofstream *out, char *path);
    … …  
    286288  void OutputListOfBonds(ofstream *out) const;
    287289  bool OutputXYZ(ofstream *out) const;
     290  bool OutputTrajectoriesXYZ(ofstream *out);
    288291  bool Checkout(ofstream *out) const;
    289292
    … …  
    354357    char *configname;
    355358    bool FastParsing;
     359    double Deltat;
    356360   
    357361    private:
    … …  
    375379   
    376380    int MaxOuterStep;
    377     double Deltat;
    378381    int OutVisStep;
    379382    int OutSrcStep;
Note: See TracChangeset for help on using the changeset viewer.