source: src/atom_bondedparticle.cpp@ d557374

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * 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 "atom.hpp"
#include "atom_bondedparticle.hpp"
#include "bond.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "element.hpp"
#include "lists.hpp"

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

/** Destructor of class BondedParticle.
 */
BondedParticle::~BondedParticle()
{
  BondList::iterator Runner;
  for (std::vector<BondList>::iterator iter = ListOfBonds.begin();
      !ListOfBonds.empty();
      iter = ListOfBonds.begin()) {
    while (!(*iter).empty()) {
      Runner = (*iter).begin();
      removewithoutcheck(*Runner);
    }
    ListOfBonds.erase(iter);
  }
};

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

/** Prints all bonds of this atom with total degree.
 */
void BondedParticle::OutputBondOfAtom() const
{
  const BondList& ListOfBonds = getListOfBonds();
  DoLog(4) && (Log() << Verbose(4) << "Atom " << getName() << "/" << nr << " with " << ListOfBonds.size() << " bonds: " << endl);
  int TotalDegree = 0;
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
    DoLog(4) && (Log() << Verbose(4) << **Runner << endl);
    TotalDegree += (*Runner)->BondDegree;
  }
  DoLog(4) && (Log() << Verbose(4) << " -- TotalDegree: " << TotalDegree << endl);
};

/** Output of atom::nr along with all bond partners.
 * \param *AdjacencyFile output stream
 */
void BondedParticle::OutputAdjacency(ofstream * const AdjacencyFile) const
{
  const BondList& ListOfBonds = getListOfBonds();
  *AdjacencyFile << nr << "\t";
  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
  *AdjacencyFile << endl;
};

/** 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 (nr < (*Runner)->GetOtherAtom(this)->nr)
      *BondFile << nr << "\t" << (*Runner)->GetOtherAtom(this)->nr << "\n";
};

/**
 * Adds a bond between this bonded particle and another. Does nothing if this
 * bond already exists.
 *
 * \param bonding partner
 */
void BondedParticle::addBond(BondedParticle* Partner) {
  if (IsBondedTo(Partner)) {
    return;
  }

  bond* newBond = new bond((atom*) this, (atom*) Partner, 1, 0);
  RegisterBond(newBond);
  Partner->RegisterBond(newBond);
}

/** Puts a given bond into atom::ListOfBonds.
 * \param *Binder bond to insert
 */
bool BondedParticle::RegisterBond(bond *Binder)
{
  bool status = false;
  if (Binder != NULL) {
    if (Binder->Contains(this)) {
      //LOG(3,"INFO: Registering bond "<< *Binder << " with atom " << *this << " at step " << _step);
      BondList& ListOfBonds = getListOfBonds();
      ListOfBonds.push_back(Binder);
      status = true;
    } else {
      DoeLog(1) && (eLog()<< Verbose(1) << *Binder << " does not contain " << *this << "." << endl);
    }
  } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Binder is " << Binder << "." << endl);
  }
  return status;
};

/** Removes a given bond from atom::ListOfBonds.
 * @param _step time step to access
 * \param *Binder bond to remove
 */
bool BondedParticle::UnregisterBond(bond *Binder)
{
  bool status = false;
  ASSERT(Binder != NULL, "BondedParticle::UnregisterBond() - Binder is NULL.");
  const int step = ContainsBondAtStep(Binder);
  if (step != -1) {
    //LOG(3,"INFO: Unregistering bond "<< *Binder << " from list " << &ListOfBonds << " of atom " << *this << " at step " << step);
    ListOfBonds[step].remove(Binder);
    status = true;
  } else {
    DoeLog(1) && (eLog()<< Verbose(1) << *Binder << " does not contain " << *this << "." << endl);
  }
  return status;
};

/** Removes all bonds from atom::ListOfBonds.
 * \note Does not do any memory de-allocation.
 */
void BondedParticle::UnregisterAllBond(const unsigned int _step)
{
  BondList& ListOfBonds = getListOfBondsAtStep(_step);
  ListOfBonds.clear();
};

/** 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 *Binder)
{
  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()
{
  int NoBonds = 0;
  int OtherNoBonds = 0;
  int FalseBondDegree = 0;
  atom *OtherWalker = NULL;
  bond *CandidateBond = NULL;

  NoBonds = CountBonds();
  //Log() << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
  if ((int)(getType()->getNoValenceOrbitals()) > NoBonds) { // 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();
      //Log() << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << OtherNoBonds << "?" << endl;
      if ((int)(OtherWalker->getType()->getNoValenceOrbitals()) > OtherNoBonds) { // 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
          CandidateBond = (*Runner);
          //Log() << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
        }
      }
    }
    if ((CandidateBond != NULL)) {
      CandidateBond->BondDegree++;
      //Log() << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
    } else {
      DoeLog(2) && (eLog()<< Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl);
      FalseBondDegree++;
    }
  }
  return FalseBondDegree;
};

/** 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)->BondDegree;
  return NoBonds;
};

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

  const BondList& ListOfBonds = getListOfBonds();
  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;
}