/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 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 .
*/
/*
* ChangeBondAngleAction.cpp
*
* Created on: Nov 14, 2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include "Actions/MoleculeAction/ChangeBondAngleAction.hpp"
#include
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "LinearAlgebra/Plane.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "molecule.hpp"
#include "World.hpp"
#include "WorldTime.hpp"
using namespace boost::assign;
using namespace MoleCuilder;
// and construct the stuff
#include "ChangeBondAngleAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */
const std::vector sortIndicesFromCentralAtom(const std::vector< atom *> &atoms)
{
std::vector indices;
{
std::vector matches;
size_t index = 0;
for (; index < 3; ++index)
if ((atoms[index]->IsBondedTo(WorldTime::getTime(), atoms[(index+1)%3]))
&& (atoms[index]->IsBondedTo(WorldTime::getTime(), atoms[(index+2)%3])))
matches.push_back(index);
if (matches.size() != 1) {
ELOG(1, "The three atoms are not linearly bonded.");
} else {
LOG(1, "INFO: Picking " << *atoms[ matches[0] ] << " as the central atom.");
indices += matches[0], (matches[0]+1)%3, (matches[0]+2)%3;
}
}
return indices;
}
void ChangeBondAngleSymmetrically(
const std::vector &indices,
const double newangle,
const std::vector &atoms)
{
const Vector firstPosition = atoms[ indices[1] ]->getPosition();
const Vector secondPosition = atoms[ indices[2] ]->getPosition();
// create the bond plane and mid-distance
const double normaldistance = firstPosition.distance(secondPosition);
const Vector NormalVector = (1. / normaldistance) * (firstPosition - secondPosition);
const Vector OffsetVector = atoms[indices[0]]->getPosition();
Plane midplane(NormalVector, OffsetVector);
// go through the molecule and stretch each atom relative two plane
for (size_t index = 1; index < 3; ++index) {
const Vector &position = atoms[indices[index]]->getPosition();
const Vector planespot = midplane.getClosestPoint(position);
const double opposite_side = position.distance(planespot);
const double adjacent_side = OffsetVector.distance(planespot);
const double hypotenuse = position.distance(OffsetVector);
if (fabs(hypotenuse) > MYEPSILON) {
const double angle = acos(adjacent_side / hypotenuse);
LOG(1, "INFO: Old angle is " << (180. / M_PI) * (2. * angle) << ".");
const double new_opposite_side = hypotenuse * sin(newangle);
const double new_adjacent_side = hypotenuse * cos(newangle);
Vector newposition = (new_opposite_side / opposite_side) * (position - planespot);
newposition += (new_adjacent_side / adjacent_side) * (planespot - OffsetVector);
newposition += OffsetVector;
atoms[indices[index]]->setPosition(newposition);
}
}
}
Action::state_ptr MoleculeChangeBondAngleAction::performCall()
{
// check precondition: three atoms
const std::vector< atom *> atoms = World::getInstance().getSelectedAtoms();
if (atoms.size() != 3) {
ELOG(1, "Exactly three atoms must be selected.");
return Action::failure;
}
// check precondition: linearly bonded atoms
const std::vector indices = sortIndicesFromCentralAtom(atoms);
if (indices.size() != 3)
return Action::failure;
const molecule *mol = atoms[0]->getMolecule();
if ((mol != atoms[1]->getMolecule()) || (mol != atoms[2]->getMolecule())) {
ELOG(1, "The two selected atoms must belong to the same molecule.");
return Action::failure;
}
// gather undo info
const Vector firstPosition = atoms[ indices[1] ]->getPosition();
const Vector secondPosition = atoms[ indices[2] ]->getPosition();
// change the bond angle
const double newangle = .5* (M_PI / 180.) * params.bondangle.get();
LOG(1, "INFO: New angle is " << (180. / M_PI) * (2. * newangle) << ".");
ChangeBondAngleSymmetrically(indices, newangle, atoms);
// create undo information
MoleculeChangeBondAngleState *UndoState = new MoleculeChangeBondAngleState(
indices[1], indices[2],
firstPosition, secondPosition,
atoms[ indices[1] ]->getPosition(), atoms[ indices[2] ]->getPosition(),
params);
return Action::state_ptr(UndoState);
}
Action::state_ptr MoleculeChangeBondAngleAction::performUndo(Action::state_ptr _state) {
MoleculeChangeBondAngleState *state = assert_cast(_state.get());
// undo both position changes
atom * const first = World::getInstance().getAtom(AtomById(state->firstId));
atom * const second = World::getInstance().getAtom(AtomById(state->secondId));
first->setPosition(state->firstOldPosition);
second->setPosition(state->secondOldPosition);
return Action::state_ptr(_state);
}
Action::state_ptr MoleculeChangeBondAngleAction::performRedo(Action::state_ptr _state){
MoleculeChangeBondAngleState *state = assert_cast(_state.get());
// redo both position changes
atom * const first = World::getInstance().getAtom(AtomById(state->firstId));
atom * const second = World::getInstance().getAtom(AtomById(state->secondId));
first->setPosition(state->firstNewPosition);
second->setPosition(state->secondNewPosition);
return Action::state_ptr(_state);
}
bool MoleculeChangeBondAngleAction::canUndo() {
return true;
}
bool MoleculeChangeBondAngleAction::shouldUndo() {
return true;
}
/** =========== end of function ====================== */