source: src/Atom/atom_bondedparticle.cpp@ 2d5137

/*
 * 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/>.
 */

/*
 * atom_bondedparticle.cpp
 *
 *  Created on: Oct 19, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <algorithm>
#include <boost/bind.hpp>

#include "atom.hpp"
#include "atom_bondedparticle.hpp"
#include "Bond/bond.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Element/element.hpp"
#include "WorldTime.hpp"

/** Constructor of class BondedParticle.
 */
BondedParticle::BondedParticle()
{
  ListOfBonds.push_back(BondList());
};

/** Destructor of class BondedParticle.
 */
BondedParticle::~BondedParticle()
{
  removeAllBonds();
};

/** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
 * \param *file output stream
 */
void BondedParticle::OutputOrder(ofstream *file) const
{
  *file << getNr() << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
  //LOG(2, "Storing: " << getNr() << "\t" << (int)AdaptiveOrder << "\t" << MaxOrder << ".");
};

/** Prints all bonds of this atom with total degree.
 */
void BondedParticle::OutputBondOfAtom(std::ostream &ost) const
{
  const BondList& ListOfBonds = getListOfBonds();
  ost << "Atom " << getName() << "/" << getNr() << " with " << ListOfBonds.size() << " bonds: ";
  int TotalDegree = 0;
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
    ost << **Runner << "\t";
    TotalDegree += (*Runner)->getDegree();
  }
  ost << " -- TotalDegree: " << TotalDegree;
};

/** Output of atom::Nr along each bond partner per line.
 * Only bonds are printed where atom::Nr is smaller than the one of the bond partner.
 * \param *AdjacencyFile output stream
 */
void BondedParticle::OutputBonds(ofstream * const BondFile) const
{
  const BondList& ListOfBonds = getListOfBonds();
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
    if (getNr() < (*Runner)->GetOtherAtom(this)->getNr())
      *BondFile << getNr() << "\t" << (*Runner)->GetOtherAtom(this)->getNr() << "\n";
};

/**
 * Adds a bond between this bonded particle and another. Returns present instance if this
 * bond already exists.
 *
 * @param _step time step to access
 * @param bonding partner
 * @return pointer to created bond or to already present bonds
 */
bond::ptr BondedParticle::addBond(const unsigned int _step, BondedParticle* Partner)
{
  const BondList &bondlist = getListOfBondsAtStep(_step);
  for (BondList::const_iterator runner = bondlist.begin();
      runner != bondlist.end();
      runner++) {
    if ((*runner)->Contains(Partner))
      return *runner;
  }

  bond::ptr newBond(new bond((atom*) this, (atom*) Partner, 1));
  RegisterBond(_step, newBond);
  Partner->RegisterBond(_step, newBond);

  return newBond;
}

/**
 * Adds a bond between this bonded particle and another. Returns present instance if this
 * bond already exists.
 *
 * @param bonding partner
 * @return pointer to created bond or to already present bonds
 */
bond::ptr BondedParticle::addBond(BondedParticle* Partner)
{
  return addBond(WorldTime::getTime(), Partner);
}

/** Helper function to find the time step to a given bond in \a Binder.
 *
 * \param Binder bond to look for
 * \return ListOfBonds::size() - not found, else - step containing \a Binder
 */
unsigned int BondedParticle::findBondsStep(bond::ptr const Binder) const
{

  size_t _step = 0;
  for (;_step < ListOfBonds.size();++_step) {
    const BondList& ListOfBonds = getListOfBondsAtStep(_step);
    if (std::find(ListOfBonds.begin(), ListOfBonds.end(), Binder) != ListOfBonds.end())
      break;
  }
  return _step;
}

/** Helper function to find the iterator to a bond at a given time \a step to
 * a given bond partner in \a Partner.
 *
 * \param _step time step to look at
 * \param Partner bond partner to look for
 * \return ListOfBonds::end() - not found, else - iterator pointing \a Binder
 */
BondList::const_iterator BondedParticle::findBondPartnerAtStep(
    const unsigned int _step,
    BondedParticle * const Partner) const
{
  const BondList& ListOfBonds = getListOfBondsAtStep(_step);
  BondList::const_iterator iter = std::find_if(ListOfBonds.begin(), ListOfBonds.end(),
      boost::bind(
          static_cast<bool (bond::*)(const ParticleInfo * const) const>(&bond::Contains),
          _1,
          boost::cref(Partner)));
  return iter;
}

/** Removes a bond of this atom to a given \a Partner.
 *
 * @param _step time step
 * @param Partner bond partner
 */
void BondedParticle::removeBond(const unsigned int _step, BondedParticle * const Partner)
{
  BondList::const_iterator iter = findBondPartnerAtStep(_step, Partner);
  if (iter != getListOfBondsAtStep(_step).end()) {
    // iter becomes invalid upon first unregister,
    // hence store the bond someplace else first
    bond::ptr const Binder = *iter;
    UnregisterBond(_step, Binder);
    Partner->UnregisterBond(_step, Binder);
  } else
    ELOG(1, "BondedParticle::removeBond() - I cannot find the bond in between "
        +toString(getName())+" and "+toString(Partner->getName())+".");
}

/** Removes a bond of this atom to a given \a Partner.
 *
 * @param Partner bond partner
 */
void BondedParticle::removeBond(BondedParticle * const Partner)
{
  removeBond(WorldTime::getTime(), Partner);
}

/** Removes a bond for this atom.
 *
 * @param Binder bond to remove
 */
void BondedParticle::removeBond(bond::ptr &binder)
{
  if (binder != NULL) {
    atom * const Other = binder->GetOtherAtom(this);
    ASSERT( Other != NULL,
        "BondedParticle::removeBonds() - cannot find bond partner for "
        +toString(*binder)+".");
    // find bond at step
    unsigned int step = findBondsStep(binder);
    if (step != ListOfBonds.size()) {
      UnregisterBond(step, binder);
      Other->UnregisterBond(step, binder);
      binder.reset();
    }
  }
}

/** Removes all bonds in all timesteps and their instances, too.
 *
 */
void BondedParticle::removeAllBonds()
{
  for (size_t index = 0; index < ListOfBonds.size(); ++index)
    removeAllBonds(index);
}

/** Removes all bonds for a given \a _step and their instances, too.
 *
 * @param _step time step to access
 */
void BondedParticle::removeAllBonds(const unsigned int _step)
{
  //LOG(3,"INFO: Clearing all bonds of " << *this << ": " << ListOfBonds[_step]);
  for (BondList::iterator iter = (ListOfBonds[_step]).begin();
      !(ListOfBonds[_step]).empty();
      iter = (ListOfBonds[_step]).begin()) {
    //LOG(3,"INFO: Clearing bond (" << *iter << ") " << *(*iter) << " of list " << &ListOfBonds);
    atom * const Other = (*iter)->GetOtherAtom(this);
    ASSERT( Other != NULL,
        "BondedParticle::removeAllBonds() - cannot find bond partner for "
        +toString(**iter)+".");
    Other->UnregisterBond(_step, *iter);
    UnregisterBond(_step, *iter);
  }
}

/** Puts a given bond into atom::ListOfBonds.
 * @param _step time step to access
 * \param *Binder bond to insert
 */
bool BondedParticle::RegisterBond(const unsigned int _step, bond::ptr const Binder)
{
  bool status = false;
  if (Binder != NULL) {
    OBSERVE;
    if (Binder->Contains(this)) {
      //LOG(3,"INFO: Registering bond "<< *Binder << " with atom " << *this << " at step " << _step);
      if (ListOfBonds.size() <= _step)
        ListOfBonds.resize(_step+1);
      ListOfBonds[_step].push_back(Binder);
      if (WorldTime::getTime() == _step)
        NOTIFY(AtomObservable::BondsAdded);
      status = true;
    } else {
      ELOG(1, *Binder << " does not contain " << *this << ".");
    }
  } else {
    ELOG(1, "Binder is " << Binder << ".");
  }
  return status;
};

/** Removes a given bond from atom::ListOfBonds.
 *
 * \warning This only removes this atom not its bond partner, i.e.
 * both atoms need to call this function to fully empty a bond.
 *
 * @param _step time step to access
 * \param *Binder bond to remove
 */
bool BondedParticle::UnregisterBond(const unsigned int _step, bond::ptr const Binder)
{
  bool status = false;
  if (Binder != NULL) {
    if (Binder->Contains(this)) {
      OBSERVE;
      //LOG(0,"INFO: Unregistering bond "<< *Binder << " from list " << &ListOfBonds << " of atom " << *this << " at step " << step);
#ifndef NDEBUG
      BondList::const_iterator iter =
          std::find(ListOfBonds[_step].begin(), ListOfBonds[_step].end(), Binder);
      ASSERT( iter != ListOfBonds[_step].end(),
          "BondedParticle::UnregisterBond() - "+toString(*Binder)+" not contained at "
          +toString(_step));
#endif
      Binder->removeAtom(this);
      ListOfBonds[_step].remove(Binder);
      if (WorldTime::getTime() == _step)
        NOTIFY(AtomObservable::BondsRemoved);
      status = true;
    } else {
      ELOG(1, *Binder << " does not contain " << *this << ".");
    }
  } else {
    ELOG(1, "Binder is " << Binder << ".");
  }
  return status;
};

/** Removes all bonds of given \a _step with freeing memory.
 *
 * @param _step time step whose bonds to free
 */
void BondedParticle::ClearBondsAtStep(const unsigned int _step)
{
  removeAllBonds(_step);
}

/** Searches for the time step where the given bond \a *Binder is a bond of this particle.
 *
 * @param Binder bond to check
 * @return >=0 - first time step where bond appears, -1 - bond not present in lists
 */
int BondedParticle::ContainsBondAtStep(bond::ptr Binder) const
{
  int step = -1;
  int tempstep = 0;
  for(std::vector<BondList>::const_iterator iter = ListOfBonds.begin();
      iter != ListOfBonds.end();
      ++iter,++tempstep) {
    for (BondList::const_iterator bonditer = iter->begin();
        bonditer != iter->end();
        ++bonditer) {
      if ((*bonditer) == Binder) {
        step = tempstep;
        break;
      }
    }
    if (step != -1)
      break;
  }

  return step;
}

/** Corrects the bond degree by one at most if necessary.
 * \return number of corrections done
 */
int BondedParticle::CorrectBondDegree()
{
  OBSERVE;
  NOTIFY(AtomObservable::BondDegreeChanged);
  int NoBonds = 0;
  int OtherNoBonds = 0;
  int FalseBondDegree = 0;
  int CandidateDeltaNoBonds = -1;
  atom *OtherWalker = NULL;
  bond::ptr CandidateBond;

  NoBonds = CountBonds();
  LOG(3, "Walker " << *this << ": " << (int)this->getType()->getNoValenceOrbitals() << " > " << NoBonds << "?");
  const int DeltaNoBonds = (int)(getType()->getNoValenceOrbitals()) - NoBonds;
  if (DeltaNoBonds > 0) { // we have a mismatch, check all bonding partners for mismatch
    const BondList& ListOfBonds = getListOfBonds();
    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
      OtherWalker = (*Runner)->GetOtherAtom(this);
      OtherNoBonds = OtherWalker->CountBonds();
      const int OtherDeltaNoBonds = (int)(OtherWalker->getType()->getNoValenceOrbitals()) - OtherNoBonds;
      LOG(3, "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->getType()->getNoValenceOrbitals() << " > " << OtherNoBonds << "?");
      if (OtherDeltaNoBonds > 0) { // check if possible candidate
        const BondList& OtherListOfBonds = OtherWalker->getListOfBonds();
        if ((CandidateBond == NULL)
            || (ListOfBonds.size() > OtherListOfBonds.size()) // pick the one with fewer number of bonds first
            || ((CandidateDeltaNoBonds < 0) || (CandidateDeltaNoBonds > OtherDeltaNoBonds)) ) // pick the one closer to fulfilling its valence first
        {
          CandidateDeltaNoBonds = OtherDeltaNoBonds;
          CandidateBond = (*Runner);
          LOG(3, "New candidate is " << *CandidateBond << ".");
        }
      }
    }
    if ((CandidateBond != NULL)) {
      CandidateBond->setDegree(CandidateBond->getDegree()+1);
      LOG(2, "Increased bond degree for bond " << *CandidateBond << ".");
    } else {
      ELOG(2, "Could not find correct degree for atom " << *this << ".");
      FalseBondDegree++;
    }
  }
  return FalseBondDegree;
};

/** Sets the weight of all connected bonds to one.
 */
void BondedParticle::resetBondDegree()
{
  OBSERVE;
  NOTIFY(BondedParticle::BondDegreeChanged);
  for (std::vector<BondList>::iterator Runner = ListOfBonds.begin();
      Runner != ListOfBonds.end();
      ++Runner)
    for (BondList::iterator BondRunner = (*Runner).begin();
        BondRunner != (*Runner).end();
        ++BondRunner)
      (*BondRunner)->setDegree(1);
};

/** Counts the number of bonds weighted by bond::BondDegree.
   * @param _step time step to access
 * \param bonds times bond::BondDegree
 */
int BondedParticle::CountBonds() const
{
  int NoBonds = 0;
  const BondList& ListOfBonds = getListOfBonds();
  for (BondList::const_iterator Runner = ListOfBonds.begin();
      Runner != ListOfBonds.end();
      (++Runner))
    NoBonds += (*Runner)->getDegree();
  return NoBonds;
};

/** Checks whether there is a bond between \a this atom and the given \a *BondPartner.
 * @param _step time step to access
 * \param *BondPartner atom to check for
 * \return true - bond exists, false - bond does not exist
 */
bool BondedParticle::IsBondedTo(const unsigned int _step, BondedParticle * const BondPartner) const
{
  bool status = false;

  const BondList& ListOfBonds = getListOfBondsAtStep(_step);
  for (BondList::const_iterator runner = ListOfBonds.begin();
      runner != ListOfBonds.end();
      runner++) {
    status = status || ((*runner)->Contains(BondPartner));
  }
  return status;
};

std::ostream & BondedParticle::operator << (std::ostream &ost) const
{
  ParticleInfo::operator<<(ost);
  ost << "," << getPosition();
  return ost;
}

std::ostream & operator << (std::ostream &ost, const BondedParticle &a)
{
  a.ParticleInfo::operator<<(ost);
  ost << "," << a.getPosition();
  return ost;
}