Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/analysis_bonds.cpp

    r68f03d r112b09  
    55 *      Author: heber
    66 */
     7
     8#include "Helpers/MemDebug.hpp"
    79
    810#include "analysis_bonds.hpp"
     
    2628  Mean = 0.;
    2729
    28   atom *Walker = mol->start;
    2930  int AtomCount = 0;
    30   while (Walker->next != mol->end) {
    31     Walker = Walker->next;
    32     const int count = Walker->ListOfBonds.size();
     31  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     32    const int count = (*iter)->ListOfBonds.size();
    3333    if (Max < count)
    3434      Max = count;
     
    5151 * \param &Max maximum distance on return, 0 if no bond between the two elements
    5252 */
    53 void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, element *type1, element *type2, double &Min, double &Mean, double &Max)
     53void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max)
    5454{
    5555  Min = 2e+6;
     
    5858
    5959  int AtomNo = 0;
    60   atom *Walker = mol->start;
    61   while (Walker->next != mol->end) {
    62     Walker = Walker->next;
    63     if (Walker->type == type1)
    64       for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++)
    65         if ((*BondRunner)->GetOtherAtom(Walker)->type == type2) {
     60  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     61    if ((*iter)->type == type1)
     62      for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
     63        if ((*BondRunner)->GetOtherAtom((*iter))->type == type2) {
    6664          const double distance = (*BondRunner)->GetDistanceSquared();
    6765          if (Min > distance)
     
    124122 * \param *InterfaceElement or NULL
    125123 */
    126 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL)
    127 {
    128   atom *Walker = NULL;
    129   atom *Runner = NULL;
     124int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL)
     125{
    130126  int count = 0;
    131127  int OtherHydrogens = 0;
     
    133129  bool InterfaceFlag = false;
    134130  bool OtherHydrogenFlag = true;
    135   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    136     Walker = (*MolWalker)->start;
    137     while (Walker->next != (*MolWalker)->end) {
    138       Walker = Walker->next;
    139       for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); MolRunner++) {
    140         Runner = (*MolRunner)->start;
    141         while (Runner->next != (*MolRunner)->end) {
    142           Runner = Runner->next;
    143           if ((Walker->type->Z  == 8) && (Runner->type->Z  == 8)) {
     131  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); ++MolWalker) {
     132    molecule::iterator Walker = (*MolWalker)->begin();
     133    for(;Walker!=(*MolWalker)->end();++Walker){
     134      for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); ++MolRunner) {
     135        molecule::iterator Runner = (*MolRunner)->begin();
     136        for(;Runner!=(*MolRunner)->end();++Runner){
     137          if (((*Walker)->type->Z  == 8) && ((*Runner)->type->Z  == 8)) {
    144138            // check distance
    145             const double distance = Runner->x.DistanceSquared(Walker->x);
     139            const double distance = (*Runner)->x.DistanceSquared((*Walker)->x);
    146140            if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
    147141              // on other atom(Runner) we check for bond to interface element and
     
    151145              OtherHydrogens = 0;
    152146              InterfaceFlag = (InterfaceElement == NULL);
    153               for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) {
    154                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner);
     147              for (BondList::const_iterator BondRunner = (*Runner)->ListOfBonds.begin(); BondRunner != (*Runner)->ListOfBonds.end(); BondRunner++) {
     148                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
    155149                // if hydrogen, check angle to be greater(!) than 30 degrees
    156150                if (OtherAtom->type->Z == 1) {
    157                   const double angle = CalculateAngle(&OtherAtom->x, &Runner->x, &Walker->x);
     151                  const double angle = CalculateAngle(&OtherAtom->x, &(*Runner)->x, &(*Walker)->x);
    158152                  OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
    159153                  Otherangle += angle;
     
    176170              if (InterfaceFlag && OtherHydrogenFlag) {
    177171                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
    178                 for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
    179                   atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
     172                for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {
     173                  atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
    180174                  if (OtherAtom->type->Z == 1) {
    181175                    // check angle
    182                     if (CheckHydrogenBridgeBondAngle(Walker, OtherAtom, Runner)) {
    183                       DoLog(1) && (Log() << Verbose(1) << Walker->getName() << ", " << OtherAtom->getName() << " and " << Runner->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &Walker->x, &Runner->x)*(180./M_PI) << "." << endl);
     176                    if (CheckHydrogenBridgeBondAngle(*Walker, OtherAtom, *Runner)) {
     177                      DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &(*Walker)->x, &(*Runner)->x)*(180./M_PI) << "." << endl);
    184178                      count++;
    185179                      break;
     
    205199int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second)
    206200{
    207   atom *Walker = NULL;
    208201  int count = 0;
    209202
    210203  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    211     Walker = (*MolWalker)->start;
    212     while (Walker->next != (*MolWalker)->end) {
    213       Walker = Walker->next;
    214       if ((Walker->type == first) || (Walker->type == second)) {  // first element matches
    215         for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
    216           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
    217           if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (Walker->nr < OtherAtom->nr)) {
     204    molecule::iterator Walker = (*MolWalker)->begin();
     205    for(;Walker!=(*MolWalker)->end();++Walker){
     206      atom * theAtom = *Walker;
     207      if ((theAtom->type == first) || (theAtom->type == second)) {  // first element matches
     208        for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
     209          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
     210          if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (theAtom->nr < OtherAtom->nr)) {
    218211            count++;
    219212            DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl);
     
    240233  bool MatchFlag[2];
    241234  bool result = false;
    242   atom *Walker = NULL;
    243235  const element * ElementArray[2];
    244236  ElementArray[0] = first;
     
    246238
    247239  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    248     Walker = (*MolWalker)->start;
    249     while (Walker->next != (*MolWalker)->end) {
    250       Walker = Walker->next;
    251       if (Walker->type == second) {  // first element matches
     240    molecule::iterator Walker = (*MolWalker)->begin();
     241    for(;Walker!=(*MolWalker)->end();++Walker){
     242      atom *theAtom = *Walker;
     243      if (theAtom->type == second) {  // first element matches
    252244        for (int i=0;i<2;i++)
    253245          MatchFlag[i] = false;
    254         for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
    255           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
     246        for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
     247          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
    256248          for (int i=0;i<2;i++)
    257249            if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) {
Note: See TracChangeset for help on using the changeset viewer.