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