source: src/Bond/StretchBond.cpp@ ea30e6

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2020 Frederik Heber. 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/>.
 */

/*
 * StretchBond.cpp
 *
 *  Created on: Oct 4, 2020
 *      Author: heber
 */


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

//#include "CodePatterns/MemDebug.hpp"

#include "StretchBond.hpp"

#include <boost/bind.hpp>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"

#include "LinearAlgebra/Plane.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Descriptors/MoleculeIdDescriptor.hpp"
#include "Element/element.hpp"
#include "Graph/BondGraph.hpp"
#include "Graph/BoostGraphCreator.hpp"
#include "Graph/BreadthFirstSearchGatherer.hpp"
#include "molecule.hpp"
#include "World.hpp"

StretchBondUtil::StretchBondUtil(
    const std::vector< atom *> _atoms
    ) :
    atoms(_atoms),
    mol(NULL),
    Shift(2, zeroVec),
    atomids(2, -1),
    domain(World::getInstance().getDomain()),
    bondside_sets(2, BoostGraphHelpers::Nodeset_t())
{
  // check preconditions
  ASSERT( atoms.size() == (size_t)2,
      "StretchBondUtil::StretchBondUtil() - exactly two atoms must be selected.");

  olddistance = atoms[0]->getPosition().distance(atoms[1]->getPosition());
  LOG(1, "INFO: Old bond distance is " << olddistance << ".");

  // gather sorted ids
  atomids[0] = atoms[0]->getId();
  atomids[1] = atoms[1]->getId();
  std::sort(atomids.begin(), atomids.end());
  LOG(1, "DEBUG: Selected nodes are " << atomids);

  mol = World::getInstance().getMolecule(MoleculeById(atoms[0]->getMolecule()->getId()));
}

bool StretchBondUtil::addEdgePredicate(
    const bond &_bond,
    const std::vector<atomId_t> &_atomids)
{
  ASSERT(_atomids.size() == (size_t)2,
      "addEdgePredicate() - atomids must contain exactly two ids.");
  // do not add selected edge
  return ((_bond.leftatom->getId() != _atomids[0])
      || (_bond.rightatom->getId() != _atomids[1]));
}

bool StretchBondUtil::operator()(double newdistance)
{
  if (newdistance <= 0.) {
    // get typical distance from table
    BondGraph *&BG = World::getInstance().getBondGraph();
    newdistance = BG->GetBondLength(
        atoms[0]->getElementNo(),
        atoms[1]->getElementNo());
    LOG(1, "INFO: Using typical bond distance of " << newdistance
        << " between " << atoms[0]->getElement().getSymbol() << " and "
        << atoms[1]->getElement().getSymbol() << ".");
  }

  const Vector NormalVector = (atoms[0]->getPosition() - atoms[1]->getPosition())* (1./olddistance);
  const double shift = 0.5*(newdistance - olddistance);
  Shift[0] = shift * NormalVector;
  Shift[1] = -shift * NormalVector;

 if (mol != atoms[1]->getMolecule()) {
    ELOG(1, "The two selected atoms must belong to the same molecule.");
    return false;
  }

  // Assume the selected bond splits the molecule into two parts, each one on
  // either side of the bond. We need to perform a BFS from each bond partner
  // not using the selected bond. Therefrom, we obtain two sets of atoms/nodes.
  // If both are disjoint, the bond is not contained in a cycle and we simply
  // shift either set as desired. If not, then we simply shift each atom,
  // leaving the other positions untouched.

  // get nodes on either side of selected bond via BFS discovery
  BoostGraphCreator BGcreator;
  BGcreator.createFromMolecule(*mol,
      boost::bind(addEdgePredicate, _1, boost::ref(atomids)));
  BreadthFirstSearchGatherer NodeGatherer(BGcreator);
  for(size_t j=0;j<2;++j) {
    bondside_sets[j] = NodeGatherer(atoms[j]->getId());
    std::sort(bondside_sets[j].begin(), bondside_sets[j].end());
  }

  // simple test whether bond has split the system in two disjoint sets or not
  bool isCyclic = false;
  if ((bondside_sets[0].size() + bondside_sets[1].size()) > BGcreator.getNumVertices()) {
      // Check whether there are common nodes in each set of distances
    if (BoostGraphHelpers::isCommonNodeInVector(bondside_sets[0], bondside_sets[1])) {
      ELOG(2, "Sets contain common node, hence bond must have been by cyclic."
          << " Shifting only bond partners.");
      for(size_t j=0;j<2;++j) {
        bondside_sets[j].clear();
        bondside_sets[j].push_back( atomids[j] );
        const Vector &position = atoms[j]->getPosition();
        atoms[j]->setPosition( domain.enforceBoundaryConditions(position+Shift[j]) );
      }
      isCyclic = true;
    }
  }

  // go through the molecule and stretch each atom in either set of nodes
  if (!isCyclic) {
    for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
      const Vector &position = (*iter)->getPosition();
      // for each atom determine in which set of nodes it is and shift accordingly
      const atomId_t &atomid = (*iter)->getId();
      if (std::binary_search(bondside_sets[0].begin(), bondside_sets[0].end(), atomid)) {
        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[0]) );
      } else if (std::binary_search(bondside_sets[1].begin(), bondside_sets[1].end(), atomid)) {
        (*iter)->setPosition( domain.enforceBoundaryConditions(position+Shift[1]) );
      } else {
        ELOG(1, "Atom " << *iter << " is not contained on either side of bond? Undoing done shifts");
        // Have to undo shifts
        for (size_t i=0;i<2;++i) {
          for (BoostGraphHelpers::Nodeset_t::const_iterator iter = bondside_sets[i].begin();
              iter != bondside_sets[i].end(); ++iter) {
            atom &walker = *World::getInstance().getAtom(AtomById(*iter));
            const Vector &position = walker.getPosition();
            walker.setPosition( domain.enforceBoundaryConditions(position-Shift[i]) );
          }
        }
        return false;
      }
    }
  }

  return true;
}