source: src/Analysis/analysis_bonds.cpp@ a58c16

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * analysis_bonds.cpp
 *
 *  Created on: Nov 7, 2009
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include "analysis_bonds.hpp"
#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Element/element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Verbose.hpp"
#include "CodePatterns/Log.hpp"
#include "molecule.hpp"
#include "MoleculeListClass.hpp"

/** Calculates the min, mean and maximum bond counts for the given molecule.
 * \param *mol molecule with atoms and atom::ListOfBonds
 * \param &Min minimum count on return
 * \param &Mean mean count on return
 * \param &Max maximum count on return
 */
void GetMaxMinMeanBondCount(const molecule * const mol, double &Min, double &Mean, double &Max)
{
  Min = 2e+6;
  Max = -2e+5;
  Mean = 0.;

  int AtomCount = 0;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    const BondList& ListOfBonds = (*iter)->getListOfBonds();
    const int count = ListOfBonds.size();
    if (Max < count)
      Max = count;
    if (Min > count)
      Min = count;
    Mean += count;
    AtomCount++;
  }
  if (((int)Mean % 2) != 0)
    ELOG(1, "Something is wrong with the bond structure, the number of bonds is not even!");
  Mean /= (double)AtomCount;
};

/** Calculates the min and max bond distance of all atoms of two given elements.
 * \param *mol molecule with atoms
 * \param *type1 one element
 * \param *type2 other element
 * \param &Min minimum distance on return, 0 if no bond between the two elements
 * \param &Mean mean distance (i.e. sum of distance for matching element pairs, divided by number) on return, 0 if no bond between the two elements
 * \param &Max maximum distance on return, 0 if no bond between the two elements
 */
void MinMeanMaxBondDistanceBetweenElements(const molecule *mol, const element *type1, const element *type2, double &Min, double &Mean, double &Max)
{
  Min = 2e+6;
  Mean = 0.;
  Max = -2e+6;

  int AtomNo = 0;
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    if ((*iter)->getType() == type1) {
      const BondList& ListOfBonds = (*iter)->getListOfBonds();
      for (BondList::const_iterator BondRunner = ListOfBonds.begin();
          BondRunner != ListOfBonds.end();
          BondRunner++)
        if ((*BondRunner)->GetOtherAtom((*iter))->getType() == type2) {
          const double distance = (*BondRunner)->GetDistanceSquared();
          if (Min > distance)
            Min = distance;
          if (Max < distance)
            Max = distance;
          Mean += sqrt(distance);
          AtomNo++;
        }
    }
  }
  if (Max < 0) {
    Max = Min = 0.;
  } else {
    Max = sqrt(Max);
    Min = sqrt(Min);
    Mean = Mean/(double)AtomNo;
  }
};

/** Calculate the angle between \a *first and \a *origin and \a *second and \a *origin.
 * \param *first first Vector
 * \param *origin origin of angle taking
 * \param *second second Vector
 * \return angle between \a *first and \a *second, both relative to origin at \a *origin.
 */
double CalculateAngle(const Vector &first, const Vector &central, const Vector &second)
{
  Vector OHBond;
  Vector OOBond;

  OHBond = first - central;
  OOBond = second - central;
  const double angle = OHBond.Angle(OOBond);
  return angle;
};

/** Checks whether the angle between \a *Oxygen and \a *Hydrogen and \a *Oxygen and \a *OtherOxygen is less than 30 degrees.
 * Note that distance criterion is not checked.
 * \param *Oxygen first oxygen atom, bonded to \a *Hydrogen
 * \param *Hydrogen hydrogen bonded to \a *Oxygen
 * \param *OtherOxygen other oxygen atom
 * \return true - angle criteria fulfilled, false - criteria not fulfilled, angle greater than 30 degrees.
 */
bool CheckHydrogenBridgeBondAngle(const atom & Oxygen, const atom & Hydrogen, const atom & OtherOxygen)
{
  Info FunctionInfo(__func__);

  // check angle
  const double angle = CalculateAngle(
      Hydrogen.getPosition(),
      Oxygen.getPosition(),
      OtherOxygen.getPosition());
  LOG(3, "INFO: Hydrogen bridge bond angle is " << angle << ", < " << M_PI*(30./180.) << "?");
  if (angle < M_PI*(30./180.)) {
    return true;
  } else {
    return false;
  }
};

/** Counts the number of hydrogen bridge bonds.
 * With \a *InterfaceElement an extra element can be specified that identifies some boundary.
 * Then, counting is for the h-bridges that connect to interface only.
 * \param *molecules molecules to count bonds
 * \param *InterfaceElement or NULL
 * \param *Interface2Element or NULL
 */
int CountHydrogenBridgeBonds(MoleculeListClass *molecules, const element * InterfaceElement = NULL, const element * Interface2Element = NULL)
{
  Info FunctionInfo(__func__);

  int count = 0;
  int OtherHydrogens = 0;
  double Otherangle = 0.;
  bool InterfaceFlag = false;
  bool Interface2Flag = false;
  bool OtherHydrogenFlag = true;
  LinkedCell::LinkedCell_View LC = World::getInstance().getLinkedCell(HBRIDGEDISTANCE);

  // go through every molecule
  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();
      MolWalker != molecules->ListOfMolecules.end();
      ++MolWalker) {
    LOG(2, "INFO: Current molecule is " << (*MolWalker)->getName() << ".");

    // go through every atom
    typedef std::set<const molecule *> Moleculeset;
    for(molecule::const_iterator Walker = const_cast<const molecule *>(*MolWalker)->begin();
        Walker != const_cast<const molecule *>(*MolWalker)->end();
        ++Walker) {
      // go through every oxygen
      if ((*Walker)->getType()->getAtomicNumber() == 8) {
        LOG(2, "INFO: Current oxygen atom is " << *(*Walker) << ".");

        // get all its neighbors
        LinkedCell::LinkedList NeighborList = LC.getAllNeighbors(HBRIDGEDISTANCE, (*Walker)->getPosition());
        // go through each candidate and gather the molecules of all other oxygens
        Moleculeset MoleculeNeighbors;
        for(LinkedCell::LinkedList::const_iterator Runner = NeighborList.begin();
            Runner != NeighborList.end(); ++Runner) {
          const atom * const OtherAtom = dynamic_cast<const atom *>(*Runner);
          if ((OtherAtom->getType()->getAtomicNumber() == 8) &&
              (OtherAtom->getMolecule() != (*MolWalker))) {
            LOG(3, "INFO: Possible neighboring molecule is " << OtherAtom->getMolecule()->getName() << ".");
            MoleculeNeighbors.insert(OtherAtom->getMolecule());
          }
        }

        // now go through the molecules
        for (Moleculeset::const_iterator moliter = MoleculeNeighbors.begin();
            moliter != MoleculeNeighbors.end();
            ++moliter) {
          LOG(2, "INFO: Current other molecule is " << (*moliter)->getName() << ".");

          // go through every other atom
          for(molecule::const_iterator Runner = (*moliter)->begin();
              Runner != (*moliter)->end();
              ++Runner) {
            // go through each oxygen
            if ((*Runner)->getType()->getAtomicNumber() == 8) {

              // check distance
              const double distance = (*Runner)->DistanceSquared(*(*Walker));
              if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) {
                LOG(2, "INFO: Distance between oxygen atom "
                    << (*Walker)->getName() << " and  "
                    << (*Runner)->getName() << " is "
                    << sqrt(distance) << ".");
                // distance >0 means  different atoms
                // on other atom(Runner) we check for bond to interface element and
                // check that O-O line is not in between the shanks of the two connected hydrogens (Otherangle > 104.5)
                OtherHydrogenFlag = true;
                Otherangle = 0.;
                OtherHydrogens = 0;
                InterfaceFlag = (InterfaceElement == NULL);
                Interface2Flag = (Interface2Element == NULL);
                const BondList& ListOfBonds = (*Runner)->getListOfBonds();
                for (BondList::const_iterator BondRunner = ListOfBonds.begin();
                    BondRunner != ListOfBonds.end();
                    BondRunner++) {
                  atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
                  // if hydrogen, check angle to be greater(!) than 30 degrees
                  if (OtherAtom->getType()->getAtomicNumber() == 1) {
                    const double angle = CalculateAngle(OtherAtom->getPosition(), (*Runner)->getPosition(), (*Walker)->getPosition());
                    OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
                    Otherangle += angle;
                    OtherHydrogens++;
                  }
                  InterfaceFlag = InterfaceFlag || (OtherAtom->getType() == InterfaceElement);
                  Interface2Flag = Interface2Flag || (OtherAtom->getType() == Interface2Element);
                }
                LOG(1, "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens.");
                switch (OtherHydrogens) {
                  case 0:
                  case 1:
                    break;
                  case 2:
                    OtherHydrogenFlag = OtherHydrogenFlag && (Otherangle > M_PI*(104.5/180.) + MYEPSILON);
                    break;
                  default: // 3 or more hydrogens ...
                    OtherHydrogenFlag = false;
                    break;
                }
                if (InterfaceFlag && Interface2Flag && OtherHydrogenFlag) {
                  // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
                  const BondList& ListOfBonds = (*Walker)->getListOfBonds();
                  for (BondList::const_iterator BondRunner = ListOfBonds.begin();
                      BondRunner != ListOfBonds.end();
                      BondRunner++) {
                    atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
                    if (OtherAtom->getType()->getAtomicNumber() == 1) {
                      // check angle
                      if (CheckHydrogenBridgeBondAngle(*(*Walker), *OtherAtom, *(*Runner))) {
                        count++;
                        break;
                      }
                    }
                  }
                }
              }
            }
          } // end go through molecules
        } // end gather molecules
      } // end go through every oxygen
    } // end go through every atom
  }
  return count;
}

/** Counts the number of bonds between two given elements.
 * \param *molecules list of molecules with all atoms
 * \param *first pointer to first element
 * \param *second pointer to second element
 * \return number of found bonds (\a *first-\a *second)
 */
int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second)
{
  int count = 0;

  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    molecule::iterator Walker = (*MolWalker)->begin();
    for(;Walker!=(*MolWalker)->end();++Walker){
      atom * theAtom = *Walker;
      if ((theAtom->getType() == first) || (theAtom->getType() == second)) {  // first element matches
        const BondList& ListOfBonds = theAtom->getListOfBonds();
        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
            BondRunner != ListOfBonds.end();
            BondRunner++) {
          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
          if (((OtherAtom->getType() == first) || (OtherAtom->getType() == second)) && (theAtom->getNr() < OtherAtom->getNr())) {
            count++;
            LOG(1, *first << "-" << *second << " bond found between " << *Walker << " and " << *OtherAtom << ".");
          }
        }
      }
    }
  }
  return count;
};

/** Counts the number of bonds between three given elements.
 * Note that we do not look for arbitrary sequence of given bonds, but \a *second will be the central atom and we check
 * whether it has bonds to both \a *first and \a *third.
 * \param *molecules list of molecules with all atoms
 * \param *first pointer to first element
 * \param *second pointer to second element
 * \param *third pointer to third element
 * \return number of found bonds (\a *first-\a *second-\a *third, \a *third-\a *second-\a *first, respectively)
 */
int CountBondsOfThree(MoleculeListClass * const molecules, const element * const first, const element * const second, const element * const third)
{
  int count = 0;
  bool MatchFlag[2];
  bool result = false;
  const element * ElementArray[2];
  ElementArray[0] = first;
  ElementArray[1] = third;

  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    molecule::iterator Walker = (*MolWalker)->begin();
    for(;Walker!=(*MolWalker)->end();++Walker){
      atom *theAtom = *Walker;
      if (theAtom->getType() == second) {  // first element matches
        for (int i=0;i<2;i++)
          MatchFlag[i] = false;
        const BondList& ListOfBonds = theAtom->getListOfBonds();
        for (BondList::const_iterator BondRunner = ListOfBonds.begin();
            BondRunner != ListOfBonds.end();
            BondRunner++) {
          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
          for (int i=0;i<2;i++)
            if ((!MatchFlag[i]) && (OtherAtom->getType() == ElementArray[i])) {
              MatchFlag[i] = true;
              break;  // each bonding atom can match at most one element we are looking for
            }
        }
        result = true;
        for (int i=0;i<2;i++) // gather results
          result = result && MatchFlag[i];
        if (result) { // check results
          count++;
          LOG(1, *first << "-" << *second << "-" << *third << " bond found at " << *Walker << ".");
        }
      }
    }
  }
  return count;
};