Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/molecule.cpp

    ra67d19 re359a8  
    44 *
    55 */
    6 
    7 #include <cstring>
    86
    97#include "atom.hpp"
     
    2321#include "tesselation.hpp"
    2422#include "vector.hpp"
    25 #include "World.hpp"
    2623
    2724/************************************* Functions for class molecule *********************************/
     
    4643  for(int i=MAX_ELEMENTS;i--;)
    4744    ElementsInMolecule[i] = 0;
    48   strcpy(name,World::get()->DefaultName);
     45  cell_size[0] = cell_size[2] = cell_size[5]= 20.;
     46  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
     47  strcpy(name,"none");
    4948};
    5049
     
    158157  double *matrix = NULL;
    159158  bond *Binder = NULL;
    160   double * const cell_size = World::get()->cell_size;
    161159
    162160//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
     
    194192  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
    195193  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);
     194    eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    197195    return false;
    198196    BondRescale = bondlength;
     
    237235            SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
    238236          } else {
    239             DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);
     237            eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;
    240238          }
    241239        }
     
    274272      bondangle = TopOrigin->type->HBondAngle[1];
    275273      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);
     274        eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
    277275        return false;
    278276        bondangle = 0;
     
    396394      break;
    397395    default:
    398       DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);
     396      eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;
    399397      AllWentWell = false;
    400398      break;
     
    429427  input = new istringstream(line);
    430428  *input >> NumberOfAtoms;
    431   DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);
     429  Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
    432430  getline(xyzfile,line,'\n'); // Read comment
    433   DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);
     431  Log() << Verbose(1) << "Comment: " << line << endl;
    434432
    435433  if (MDSteps == 0) // no atoms yet present
     
    447445    Walker->type = elemente->FindElement(shorthand);
    448446    if (Walker->type == NULL) {
    449       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
     447      eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";
    450448      Walker->type = elemente->FindElement(1);
    451449    }
     
    543541    add(Binder, last);
    544542  } 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);
     543    eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;
    546544  }
    547545  return Binder;
     
    555553bool molecule::RemoveBond(bond *pointer)
    556554{
    557   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
     555  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
    558556  pointer->leftatom->RegisterBond(pointer);
    559557  pointer->rightatom->RegisterBond(pointer);
     
    569567bool molecule::RemoveBonds(atom *BondPartner)
    570568{
    571   //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);
     569  //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
    572570  BondList::const_iterator ForeRunner;
    573571  while (!BondPartner->ListOfBonds.empty()) {
     
    589587  else
    590588    molname = filename; // contains no slashes
    591   const char *endname = strchr(molname, '.');
     589  char *endname = strchr(molname, '.');
    592590  if ((endname == NULL) || (endname < molname))
    593591    length = strlen(molname);
     
    603601void molecule::SetBoxDimension(Vector *dim)
    604602{
    605   double * const cell_size = World::get()->cell_size;
    606603  cell_size[0] = dim->x[0];
    607604  cell_size[1] = 0.;
     
    622619    AtomCount--;
    623620  } 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);
     621    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    625622  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    626623    ElementCount--;
     
    640637    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
    641638  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);
     639    eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    643640  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    644641    ElementCount--;
     
    665662    return walker;
    666663  } else {
    667     DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
     664    Log() << Verbose(0) << "Atom not found in list." << endl;
    668665    return NULL;
    669666  }
     
    681678    //mol->Output((ofstream *)&cout);
    682679    //Log() << Verbose(0) << "===============================================" << endl;
    683     DoLog(0) && (Log() << Verbose(0) << text);
     680    Log() << Verbose(0) << text;
    684681    cin >> No;
    685682    ion = this->FindAtom(No);
     
    694691bool molecule::CheckBounds(const Vector *x) const
    695692{
    696   double * const cell_size = World::get()->cell_size;
    697693  bool result = true;
    698694  int j =-1;
     
    770766void molecule::OutputListOfBonds() const
    771767{
    772   DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
     768  Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;
    773769  ActOnAllAtoms (&atom::OutputBondOfAtom );
    774   DoLog(0) && (Log() << Verbose(0) << endl);
     770  Log() << Verbose(0) << endl;
    775771};
    776772
     
    829825  }
    830826  if ((AtomCount == 0) || (i != AtomCount)) {
    831     DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
     827    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
    832828    AtomCount = i;
    833829
     
    845841        Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
    846842        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);
     843        Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;
    848844        i++;
    849845      }
    850846    } else
    851       DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
     847      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
    852848  }
    853849};
     
    909905  bool result = true; // status of comparison
    910906
    911   DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
     907  Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
    912908  /// first count both their atoms and elements and update lists thereby ...
    913909  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
     
    921917  if (result) {
    922918    if (AtomCount != OtherMolecule->AtomCount) {
    923       DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
     919      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
    924920      result = false;
    925921    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
     
    928924  if (result) {
    929925    if (ElementCount != OtherMolecule->ElementCount) {
    930       DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);
     926      Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
    931927      result = false;
    932928    } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;
     
    940936    }
    941937    if (flag < MAX_ELEMENTS) {
    942       DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);
     938      Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;
    943939      result = false;
    944940    } else Log() << Verbose(4) << "ElementsInMolecule match." << endl;
     
    946942  /// then determine and compare center of gravity for each molecule ...
    947943  if (result) {
    948     DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
     944    Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;
    949945    DeterminePeriodicCenter(CenterOfGravity);
    950946    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    951     DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: ");
     947    Log() << Verbose(5) << "Center of Gravity: ";
    952948    CenterOfGravity.Output();
    953     DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ");
     949    Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
    954950    OtherCenterOfGravity.Output();
    955     DoLog(0) && (Log() << Verbose(0) << endl);
     951    Log() << Verbose(0) << endl;
    956952    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
    957       DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
     953      Log() << Verbose(4) << "Centers of gravity don't match." << endl;
    958954      result = false;
    959955    }
     
    962958  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
    963959  if (result) {
    964     DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
     960    Log() << Verbose(5) << "Calculating distances" << endl;
    965961    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
    966962    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
     
    969965
    970966    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    971     DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
     967    Log() << Verbose(5) << "Sorting distances" << endl;
    972968    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
    973969    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
     
    975971    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
    976972    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
    977     DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
     973    Log() << Verbose(5) << "Combining Permutation Maps" << endl;
    978974    for(int i=AtomCount;i--;)
    979975      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    980976
    981977    /// ... 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);
     978    Log() << Verbose(4) << "Comparing distances" << endl;
    983979    flag = 0;
    984980    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);
     981      Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
    986982      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
    987983        flag = 1;
     
    999995  }
    1000996  /// return pointer to map if all distances were below \a threshold
    1001   DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
     997  Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
    1002998  if (result) {
    1003     DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
     999    Log() << Verbose(3) << "Result: Equal." << endl;
    10041000    return PermutationMap;
    10051001  } else {
    1006     DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
     1002    Log() << Verbose(3) << "Result: Not equal." << endl;
    10071003    return NULL;
    10081004  }
     
    10191015{
    10201016  atom *Walker = NULL, *OtherWalker = NULL;
    1021   DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
     1017  Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
    10221018  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
    10231019  for (int i=AtomCount;i--;)
     
    10261022    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
    10271023      AtomicMap[i] = i;
    1028     DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
     1024    Log() << Verbose(4) << "Map is trivial." << endl;
    10291025  } else {
    1030     DoLog(4) && (Log() << Verbose(4) << "Map is ");
     1026    Log() << Verbose(4) << "Map is ";
    10311027    Walker = start;
    10321028    while (Walker->next != end) {
     
    10451041        }
    10461042      }
    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);
     1043      Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";
     1044    }
     1045    Log() << Verbose(0) << endl;
     1046  }
     1047  Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;
    10521048  return AtomicMap;
    10531049};
Note: See TracChangeset for help on using the changeset viewer.