Changeset 5f612ee for src/molecule.cpp


Ignore:
Timestamp:
Apr 27, 2010, 2:25:42 PM (15 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:
632bc3
Parents:
13d5a9 (diff), c695c9 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'Analysis_PairCorrelation' into StructureRefactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/World.cpp
molecuilder/src/World.hpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/log.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/periodentafel.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/Makefile.am
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/unittests/gslvectorunittest.cpp
molecuilder/src/unittests/logunittest.cpp
molecuilder/src/unittests/tesselation_boundarytriangleunittest.hpp
molecuilder/src/vector.cpp
molecuilder/tests/Tesselations/defs.in

Conflicts have been many and too numerous to listen here, just the few general cases

  • new molecule() replaced by World::getInstance().createMolecule()
  • new atom() replaced by World::getInstance().createAtom() where appropriate.
  • Some DoLog()s added interfered with changes to the message produced by Log() << Verbose(.) << ...
  • DoLog() has been erroneously added to TestRunner.cpp as well, there cout is appropriate
  • ...

Additionally, there was a bug in atom::clone(), sort was set to atom::nr of the atom to clone not of the clone itself. This caused a failure of the fragmentation.

This merge has been fully checked from a clean build directory with subsequent configure,make all install and make check.
It configures, compiles and runs all test cases and the test suite without errors.

Signed-off-by: Frederik Heber <heber@…>

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/molecule.cpp

    r13d5a9 r5f612ee  
    2525#include "tesselation.hpp"
    2626#include "vector.hpp"
     27#include "World.hpp"
    2728
    2829/************************************* Functions for class molecule *********************************/
     
    5051  for(int i=MAX_ELEMENTS;i--;)
    5152    ElementsInMolecule[i] = 0;
    52   cell_size[0] = cell_size[2] = cell_size[5]= 20.;
    53   cell_size[1] = cell_size[3] = cell_size[4]= 0.;
    54   strcpy(name,"none");
     53  strcpy(name,World::getInstance().getDefaultName());
    5554};
    5655
     
    215214  double *matrix = NULL;
    216215  bond *Binder = NULL;
     216  double * const cell_size = World::getInstance().getDomain();
    217217
    218218//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
     
    250250  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
    251251  if (BondRescale == -1) {
    252     eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
     252    DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
    253253    return false;
    254254    BondRescale = bondlength;
     
    293293            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
    294294          } else {
    295             eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;
     295            DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);
    296296          }
    297297        }
     
    330330      bondangle = TopOrigin->type->HBondAngle[1];
    331331      if (bondangle == -1) {
    332         eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
     332        DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
    333333        return false;
    334334        bondangle = 0;
     
    452452      break;
    453453    default:
    454       eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;
     454      DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);
    455455      AllWentWell = false;
    456456      break;
     
    487487  input = new istringstream(line);
    488488  *input >> NumberOfAtoms;
    489   Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
     489  DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);
    490490  getline(xyzfile,line,'\n'); // Read comment
    491   Log() << Verbose(1) << "Comment: " << line << endl;
     491  DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);
    492492
    493493  if (MDSteps == 0) // no atoms yet present
     
    505505    Walker->type = elemente->FindElement(shorthand);
    506506    if (Walker->type == NULL) {
    507       eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";
     507      DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
    508508      Walker->type = elemente->FindElement(1);
    509509    }
     
    532532molecule *molecule::CopyMolecule()
    533533{
    534   molecule *copy = new molecule(elemente);
     534  molecule *copy = World::getInstance().createMolecule();
    535535  atom *LeftAtom = NULL, *RightAtom = NULL;
    536536
     
    575575 */
    576576molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
    577   molecule *copy = new molecule(elemente);
     577  molecule *copy = World::getInstance().createMolecule();
    578578
    579579  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
     
    601601    add(Binder, last);
    602602  } else {
    603     eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
     603    DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl);
    604604  }
    605605  return Binder;
     
    613613bool molecule::RemoveBond(bond *pointer)
    614614{
    615   //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
     615  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
    616616  pointer->leftatom->RegisterBond(pointer);
    617617  pointer->rightatom->RegisterBond(pointer);
     
    627627bool molecule::RemoveBonds(atom *BondPartner)
    628628{
    629   //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
     629  //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
    630630  BondList::const_iterator ForeRunner;
    631631  while (!BondPartner->ListOfBonds.empty()) {
     
    661661void molecule::SetBoxDimension(Vector *dim)
    662662{
     663  double * const cell_size = World::getInstance().getDomain();
    663664  cell_size[0] = dim->x[0];
    664665  cell_size[1] = 0.;
     
    679680    AtomCount--;
    680681  } else
    681     eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
     682    DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);
    682683  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    683684    ElementCount--;
     
    697698    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
    698699  else
    699     eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
     700    DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);
    700701  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    701702    ElementCount--;
     
    722723    return walker;
    723724  } else {
    724     Log() << Verbose(0) << "Atom not found in list." << endl;
     725    DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
    725726    return NULL;
    726727  }
     
    738739    //mol->Output((ofstream *)&cout);
    739740    //Log() << Verbose(0) << "===============================================" << endl;
    740     Log() << Verbose(0) << text;
     741    DoLog(0) && (Log() << Verbose(0) << text);
    741742    cin >> No;
    742743    ion = this->FindAtom(No);
     
    751752bool molecule::CheckBounds(const Vector *x) const
    752753{
     754  double * const cell_size = World::getInstance().getDomain();
    753755  bool result = true;
    754756  int j =-1;
     
    826828void molecule::OutputListOfBonds() const
    827829{
    828   Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
     830  DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
    829831  ActOnAllAtoms (&atom::OutputBondOfAtom );
    830   Log() << Verbose(0) << endl;
     832  DoLog(0) && (Log() << Verbose(0) << endl);
    831833};
    832834
     
    885887  }
    886888  if ((AtomCount == 0) || (i != AtomCount)) {
    887     Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
     889    DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
    888890    AtomCount = i;
    889891
     
    901903        Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
    902904        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
    903         Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
     905        DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl);
    904906        i++;
    905907      }
    906908    } else
    907       Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
     909      DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
    908910  }
    909911};
     
    965967  bool result = true; // status of comparison
    966968
    967   Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
     969  DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
    968970  /// first count both their atoms and elements and update lists thereby ...
    969971  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
     
    977979  if (result) {
    978980    if (AtomCount != OtherMolecule->AtomCount) {
    979       Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
     981      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
    980982      result = false;
    981983    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
     
    984986  if (result) {
    985987    if (ElementCount != OtherMolecule->ElementCount) {
    986       Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
     988      DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);
    987989      result = false;
    988990    } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
     
    996998    }
    997999    if (flag < MAX_ELEMENTS) {
    998       Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;
     1000      DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);
    9991001      result = false;
    10001002    } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
     
    10021004  /// then determine and compare center of gravity for each molecule ...
    10031005  if (result) {
    1004     Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;
     1006    DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
    10051007    DeterminePeriodicCenter(CenterOfGravity);
    10061008    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    1007     Log() << Verbose(5) << "Center of Gravity: ";
     1009    DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: ");
    10081010    CenterOfGravity.Output();
    1009     Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
     1011    DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ");
    10101012    OtherCenterOfGravity.Output();
    1011     Log() << Verbose(0) << endl;
     1013    DoLog(0) && (Log() << Verbose(0) << endl);
    10121014    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
    1013       Log() << Verbose(4) << "Centers of gravity don't match." << endl;
     1015      DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
    10141016      result = false;
    10151017    }
     
    10181020  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
    10191021  if (result) {
    1020     Log() << Verbose(5) << "Calculating distances" << endl;
     1022    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
    10211023    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
    10221024    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
     
    10251027
    10261028    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    1027     Log() << Verbose(5) << "Sorting distances" << endl;
     1029    DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
    10281030    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
    10291031    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
     
    10311033    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
    10321034    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
    1033     Log() << Verbose(5) << "Combining Permutation Maps" << endl;
     1035    DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
    10341036    for(int i=AtomCount;i--;)
    10351037      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    10361038
    10371039    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
    1038     Log() << Verbose(4) << "Comparing distances" << endl;
     1040    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
    10391041    flag = 0;
    10401042    for (int i=0;i<AtomCount;i++) {
    1041       Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
     1043      DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl);
    10421044      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
    10431045        flag = 1;
     
    10551057  }
    10561058  /// return pointer to map if all distances were below \a threshold
    1057   Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
     1059  DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
    10581060  if (result) {
    1059     Log() << Verbose(3) << "Result: Equal." << endl;
     1061    DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
    10601062    return PermutationMap;
    10611063  } else {
    1062     Log() << Verbose(3) << "Result: Not equal." << endl;
     1064    DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
    10631065    return NULL;
    10641066  }
     
    10751077{
    10761078  atom *Walker = NULL, *OtherWalker = NULL;
    1077   Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
     1079  DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
    10781080  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
    10791081  for (int i=AtomCount;i--;)
     
    10821084    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
    10831085      AtomicMap[i] = i;
    1084     Log() << Verbose(4) << "Map is trivial." << endl;
     1086    DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
    10851087  } else {
    1086     Log() << Verbose(4) << "Map is ";
     1088    DoLog(4) && (Log() << Verbose(4) << "Map is ");
    10871089    Walker = start;
    10881090    while (Walker->next != end) {
     
    11011103        }
    11021104      }
    1103       Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
    1104     }
    1105     Log() << Verbose(0) << endl;
    1106   }
    1107   Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;
     1105      DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
     1106    }
     1107    DoLog(0) && (Log() << Verbose(0) << endl);
     1108  }
     1109  DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl);
    11081110  return AtomicMap;
    11091111};
Note: See TracChangeset for help on using the changeset viewer.