source: src/Fragmentation/Exporters/SaturatedBond.cpp@ 0983e6

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

/*
 * SaturatedBond.cpp
 *
 *  Created on: Jul 27, 2014
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "SaturatedBond.hpp"

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

#include <cmath>
#include <iostream>

#include "Atom/atom.hpp"
#include "Element/element.hpp"

SaturatedBond::SaturatedBond(
    const bond &_bond,
    const atom& _remaining) :
    saturated_bond(_bond),
    saturated_atom(_remaining)
{
  // create bond and orthogonal vectors
  const atom &other = *saturated_bond.GetOtherAtom(&saturated_atom);
  BondVector = other.getPosition() - saturated_atom.getPosition();
  BondVector.Normalize();
  vector_a.GetOneNormalVector(BondVector);
  vector_b = BondVector;
  vector_b.VectorProduct(vector_a);

  // gather data for this element from the tables
  const double HydrogenDistance =
      saturated_atom.getElement().getHBondDistance(saturated_bond.getDegree()-1);
  if( HydrogenDistance <= 0.) {
      ELOG(1, "SaturatedBond::SaturatedBond() - negative bond distance for "
          << saturated_atom.getElement().getName() << ", defaulting to 1.");
      const_cast<double &>(HydrogenDistance) = 1.;
  }
  const double HydrogenAngle =
      saturated_atom.getElement().getHBondAngle(saturated_bond.getDegree()-1);
  if ( HydrogenAngle < 0.) {
      ELOG(1, "SaturatedBond::SaturatedBond() - negative bond angle for "
          << saturated_atom.getElement().getName() << ", defaulting to 180.");
      const_cast<double &>(HydrogenAngle) = 180.;
  }
  LOG(5, "DEBUG: Hydrogen distance is " << HydrogenDistance
      << ", angle is " << HydrogenAngle);

  // create correct lengths
  const double BondVectorLength = cos((HydrogenAngle/180.*M_PI)/2.)*HydrogenDistance;
  const double OrthogonalLength = sin((HydrogenAngle/180.*M_PI)/2.)*HydrogenDistance;
  BondVector *= BondVectorLength;
  vector_a *= OrthogonalLength;
  vector_b *= OrthogonalLength;

  LOG(5, "DEBUG: BondVector is " << BondVector
      << ", vector_a is " << vector_a
      << ", and vector_b is " << vector_b);
}

std::vector<Vector> SaturatedBond::getPositions() const
{
  std::vector<Vector> positions;
  positions.reserve(saturated_bond.getDegree());

  // return positions with respect to the bond degree
  switch (saturated_bond.getDegree())
  {
    case 1:
      {
        positions.push_back(BondVector);
      }
      break;
    case 2:
      {
        const Vector position_a = cos(alpha) * vector_a + sin(alpha) * vector_b;
        const Vector position_b = cos(alpha+M_PI) * vector_a + sin(alpha+M_PI) * vector_b;
        positions.push_back(BondVector + position_a);
        positions.push_back(BondVector + position_b);
      }
      break;
    case 3:
      {
        const Vector position_a =
            cos(alpha) * vector_a + sin(alpha) * vector_b;
        const Vector position_b =
            cos(alpha+2.*M_PI/3.) * vector_a + sin(alpha+2.*M_PI/3.) * vector_b;
        const Vector position_c =
            cos(alpha+4.*M_PI/3.) * vector_a + sin(alpha+4.*M_PI/3.) * vector_b;
        positions.push_back(BondVector + position_a);
        positions.push_back(BondVector + position_b);
        positions.push_back(BondVector + position_c);
      }
      break;
    default:
      ASSERT(0, "SaturatedBond::getPositions() - we cannot deal with bond degree "
          +toString(saturated_bond.getDegree())+" of bond "+toString(saturated_bond));
  }

  return positions;
}

std::ostream &operator<<(std::ostream &ost, const SaturatedBond& _bond)
{
  ost << _bond.saturated_bond;
  return ost;
}