Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/molecule.cpp

    ra67d19 r49e1ae  
    2323#include "tesselation.hpp"
    2424#include "vector.hpp"
    25 #include "World.hpp"
    2625
    2726/************************************* Functions for class molecule *********************************/
     
    4645  for(int i=MAX_ELEMENTS;i--;)
    4746    ElementsInMolecule[i] = 0;
    48   strcpy(name,World::get()->DefaultName);
     47  cell_size[0] = cell_size[2] = cell_size[5]= 20.;
     48  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
     49  strcpy(name,"none");
    4950};
    5051
     
    158159  double *matrix = NULL;
    159160  bond *Binder = NULL;
    160   double * const cell_size = World::get()->cell_size;
    161161
    162162//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
     
    194194  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
    195195  if (BondRescale == -1) {
    196     DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
     196    eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    197197    return false;
    198198    BondRescale = bondlength;
     
    237237            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
    238238          } else {
    239             DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);
     239            eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;
    240240          }
    241241        }
     
    274274      bondangle = TopOrigin->type->HBondAngle[1];
    275275      if (bondangle == -1) {
    276         DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);
     276        eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    277277        return false;
    278278        bondangle = 0;
     
    396396      break;
    397397    default:
    398       DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);
     398      eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;
    399399      AllWentWell = false;
    400400      break;
     
    429429  input = new istringstream(line);
    430430  *input >> NumberOfAtoms;
    431   DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);
     431  Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
    432432  getline(xyzfile,line,'\n'); // Read comment
    433   DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);
     433  Log() << Verbose(1) << "Comment: " << line << endl;
    434434
    435435  if (MDSteps == 0) // no atoms yet present
     
    447447    Walker->type = elemente->FindElement(shorthand);
    448448    if (Walker->type == NULL) {
    449       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
     449      eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";
    450450      Walker->type = elemente->FindElement(1);
    451451    }
     
    543543    add(Binder, last);
    544544  } else {
    545     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);
     545    eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
    546546  }
    547547  return Binder;
     
    555555bool molecule::RemoveBond(bond *pointer)
    556556{
    557   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
     557  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
    558558  pointer->leftatom->RegisterBond(pointer);
    559559  pointer->rightatom->RegisterBond(pointer);
     
    569569bool molecule::RemoveBonds(atom *BondPartner)
    570570{
    571   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
     571  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
    572572  BondList::const_iterator ForeRunner;
    573573  while (!BondPartner->ListOfBonds.empty()) {
     
    603603void molecule::SetBoxDimension(Vector *dim)
    604604{
    605   double * const cell_size = World::get()->cell_size;
    606605  cell_size[0] = dim->x[0];
    607606  cell_size[1] = 0.;
     
    622621    AtomCount--;
    623622  } else
    624     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);
     623    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    625624  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    626625    ElementCount--;
     
    640639    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
    641640  else
    642     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);
     641    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    643642  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    644643    ElementCount--;
     
    665664    return walker;
    666665  } else {
    667     DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
     666    Log() << Verbose(0) << "Atom not found in list." << endl;
    668667    return NULL;
    669668  }
     
    681680    //mol->Output((ofstream *)&cout);
    682681    //Log() << Verbose(0) << "===============================================" << endl;
    683     DoLog(0) && (Log() << Verbose(0) << text);
     682    Log() << Verbose(0) << text;
    684683    cin >> No;
    685684    ion = this->FindAtom(No);
     
    694693bool molecule::CheckBounds(const Vector *x) const
    695694{
    696   double * const cell_size = World::get()->cell_size;
    697695  bool result = true;
    698696  int j =-1;
     
    770768void molecule::OutputListOfBonds() const
    771769{
    772   DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
     770  Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
    773771  ActOnAllAtoms (&atom::OutputBondOfAtom );
    774   DoLog(0) && (Log() << Verbose(0) << endl);
     772  Log() << Verbose(0) << endl;
    775773};
    776774
     
    829827  }
    830828  if ((AtomCount == 0) || (i != AtomCount)) {
    831     DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
     829    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
    832830    AtomCount = i;
    833831
     
    845843        Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
    846844        sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1);
    847         DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl);
     845        Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
    848846        i++;
    849847      }
    850848    } else
    851       DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
     849      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
    852850  }
    853851};
     
    909907  bool result = true; // status of comparison
    910908
    911   DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
     909  Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
    912910  /// first count both their atoms and elements and update lists thereby ...
    913911  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
     
    921919  if (result) {
    922920    if (AtomCount != OtherMolecule->AtomCount) {
    923       DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
     921      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
    924922      result = false;
    925923    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
     
    928926  if (result) {
    929927    if (ElementCount != OtherMolecule->ElementCount) {
    930       DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);
     928      Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
    931929      result = false;
    932930    } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
     
    940938    }
    941939    if (flag < MAX_ELEMENTS) {
    942       DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);
     940      Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;
    943941      result = false;
    944942    } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
     
    946944  /// then determine and compare center of gravity for each molecule ...
    947945  if (result) {
    948     DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
     946    Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;
    949947    DeterminePeriodicCenter(CenterOfGravity);
    950948    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    951     DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: ");
     949    Log() << Verbose(5) << "Center of Gravity: ";
    952950    CenterOfGravity.Output();
    953     DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ");
     951    Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
    954952    OtherCenterOfGravity.Output();
    955     DoLog(0) && (Log() << Verbose(0) << endl);
     953    Log() << Verbose(0) << endl;
    956954    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
    957       DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
     955      Log() << Verbose(4) << "Centers of gravity don't match." << endl;
    958956      result = false;
    959957    }
     
    962960  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
    963961  if (result) {
    964     DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
     962    Log() << Verbose(5) << "Calculating distances" << endl;
    965963    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
    966964    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
     
    969967
    970968    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    971     DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
     969    Log() << Verbose(5) << "Sorting distances" << endl;
    972970    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
    973971    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
     
    975973    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
    976974    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
    977     DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
     975    Log() << Verbose(5) << "Combining Permutation Maps" << endl;
    978976    for(int i=AtomCount;i--;)
    979977      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    980978
    981979    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
    982     DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
     980    Log() << Verbose(4) << "Comparing distances" << endl;
    983981    flag = 0;
    984982    for (int i=0;i<AtomCount;i++) {
    985       DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl);
     983      Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
    986984      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
    987985        flag = 1;
     
    999997  }
    1000998  /// return pointer to map if all distances were below \a threshold
    1001   DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
     999  Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
    10021000  if (result) {
    1003     DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
     1001    Log() << Verbose(3) << "Result: Equal." << endl;
    10041002    return PermutationMap;
    10051003  } else {
    1006     DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
     1004    Log() << Verbose(3) << "Result: Not equal." << endl;
    10071005    return NULL;
    10081006  }
     
    10191017{
    10201018  atom *Walker = NULL, *OtherWalker = NULL;
    1021   DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
     1019  Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
    10221020  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
    10231021  for (int i=AtomCount;i--;)
     
    10261024    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
    10271025      AtomicMap[i] = i;
    1028     DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
     1026    Log() << Verbose(4) << "Map is trivial." << endl;
    10291027  } else {
    1030     DoLog(4) && (Log() << Verbose(4) << "Map is ");
     1028    Log() << Verbose(4) << "Map is ";
    10311029    Walker = start;
    10321030    while (Walker->next != end) {
     
    10451043        }
    10461044      }
    1047       DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
    1048     }
    1049     DoLog(0) && (Log() << Verbose(0) << endl);
    1050   }
    1051   DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl);
     1045      Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
     1046    }
     1047    Log() << Verbose(0) << endl;
     1048  }
     1049  Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;
    10521050  return AtomicMap;
    10531051};
Note: See TracChangeset for help on using the changeset viewer.