source: src/molecule.cpp@ 748fc7

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 748fc7 was 06804b, checked in by Frederik Heber <heber@…>, 10 years ago

Removed molecule::AddHydrogenReplacement() and BreadthFirstSearchAdd.

  • this capability has been transfered to SaturatedFragment, and BFSAdd has not been used for a long, long time.
  • Property mode set to 100755
File size: 38.5 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/** \file molecules.cpp
25 *
26 * Functions for the class molecule.
27 *
28 */
29
30// include config.h
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include "CodePatterns/MemDebug.hpp"
36
37#include <cstring>
38#include <boost/bind.hpp>
39#include <boost/foreach.hpp>
40
41#include <gsl/gsl_inline.h>
42#include <gsl/gsl_heapsort.h>
43
44#include "molecule.hpp"
45
46#include "Atom/atom.hpp"
47#include "Bond/bond.hpp"
48#include "Box.hpp"
49#include "CodePatterns/enumeration.hpp"
50#include "CodePatterns/Log.hpp"
51#include "config.hpp"
52#include "Descriptors/AtomIdDescriptor.hpp"
53#include "Element/element.hpp"
54#include "Graph/BondGraph.hpp"
55#include "LinearAlgebra/Exceptions.hpp"
56#include "LinearAlgebra/leastsquaremin.hpp"
57#include "LinearAlgebra/Plane.hpp"
58#include "LinearAlgebra/RealSpaceMatrix.hpp"
59#include "LinearAlgebra/Vector.hpp"
60#include "LinkedCell/linkedcell.hpp"
61#include "IdPool_impl.hpp"
62#include "Shapes/BaseShapes.hpp"
63#include "Tesselation/tesselation.hpp"
64#include "World.hpp"
65#include "WorldTime.hpp"
66
67
68/************************************* Functions for class molecule *********************************/
69
70/** Constructor of class molecule.
71 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
72 */
73molecule::molecule() :
74 Observable("molecule"),
75 MDSteps(0),
76 NoNonBonds(0),
77 NoCyclicBonds(0),
78 ActiveFlag(false),
79 IndexNr(-1),
80 NoNonHydrogen(this,boost::bind(&molecule::doCountNoNonHydrogen,this),"NoNonHydrogen"),
81 BondCount(this,boost::bind(&molecule::doCountBonds,this),"BondCount"),
82 atomIdPool(1, 20, 100),
83 last_atom(0)
84{
85 // add specific channels
86 Channels *OurChannel = new Channels;
87 NotificationChannels.insert( std::make_pair( static_cast<Observable *>(this), OurChannel) );
88 for (size_t type = 0; type < (size_t)NotificationType_MAX; ++type)
89 OurChannel->addChannel(type);
90
91 strcpy(name,World::getInstance().getDefaultName().c_str());
92};
93
94molecule *NewMolecule(){
95 return new molecule();
96}
97
98/** Destructor of class molecule.
99 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
100 */
101molecule::~molecule()
102{
103 CleanupMolecule();
104};
105
106
107void DeleteMolecule(molecule *mol){
108 delete mol;
109}
110
111// getter and setter
112const std::string molecule::getName() const{
113 return std::string(name);
114}
115
116int molecule::getAtomCount() const{
117 return atomIds.size();
118}
119
120size_t molecule::getNoNonHydrogen() const{
121 return *NoNonHydrogen;
122}
123
124int molecule::getBondCount() const{
125 return *BondCount;
126}
127
128void molecule::setName(const std::string _name){
129 OBSERVE;
130 NOTIFY(MoleculeNameChanged);
131 cout << "Set name of molecule " << getId() << " to " << _name << endl;
132 strncpy(name,_name.c_str(),MAXSTRINGSIZE);
133}
134
135void molecule::InsertLocalToGlobalId(atom * const pointer)
136{
137#ifndef NDEBUG
138 std::pair< LocalToGlobalId_t::iterator, bool > inserter =
139#endif
140 LocalToGlobalId.insert( std::make_pair(pointer->getNr(), pointer) );
141 ASSERT( inserter.second,
142 "molecule::AddAtom() - local number "+toString(pointer->getNr())+" appears twice.");
143}
144
145bool molecule::changeAtomNr(int oldNr, int newNr, atom* target){
146 OBSERVE;
147 if(atomIdPool.reserveId(newNr)){
148 NOTIFY(AtomNrChanged);
149 if (oldNr != -1) // -1 is reserved and indicates no number
150 atomIdPool.releaseId(oldNr);
151 LocalToGlobalId.erase(oldNr);
152 ASSERT (target,
153 "molecule::changeAtomNr() - given target is NULL, cannot set Nr or name.");
154 target->setNr(newNr);
155 InsertLocalToGlobalId(target);
156 setAtomName(target);
157 return true;
158 } else{
159 return false;
160 }
161}
162
163bool molecule::changeId(moleculeId_t newId){
164 // first we move ourselves in the world
165 // the world lets us know if that succeeded
166 if(World::getInstance().changeMoleculeId(id,newId,this)){
167 id = newId;
168 return true;
169 }
170 else{
171 return false;
172 }
173}
174
175
176moleculeId_t molecule::getId() const {
177 return id;
178}
179
180void molecule::setId(moleculeId_t _id){
181 id =_id;
182}
183
184const Formula &molecule::getFormula() const {
185 return formula;
186}
187
188unsigned int molecule::getElementCount() const{
189 return formula.getElementCount();
190}
191
192bool molecule::hasElement(const element *element) const{
193 return formula.hasElement(element);
194}
195
196bool molecule::hasElement(atomicNumber_t Z) const{
197 return formula.hasElement(Z);
198}
199
200bool molecule::hasElement(const string &shorthand) const{
201 return formula.hasElement(shorthand);
202}
203
204/************************** Access to the List of Atoms ****************/
205
206molecule::const_iterator molecule::erase( const_iterator loc )
207{
208 OBSERVE;
209 NOTIFY(AtomRemoved);
210 const_iterator iter = loc;
211 ++iter;
212 atom * const _atom = const_cast<atom *>(*loc);
213 atomIds.erase( _atom->getId() );
214 {
215 NOTIFY(AtomNrChanged);
216 atomIdPool.releaseId(_atom->getNr());
217 LocalToGlobalId.erase(_atom->getNr());
218 _atom->setNr(-1);
219 }
220 formula-=_atom->getType();
221 _atom->removeFromMolecule();
222 return iter;
223}
224
225molecule::const_iterator molecule::erase( atom * key )
226{
227 OBSERVE;
228 NOTIFY(AtomRemoved);
229 const_iterator iter = find(key);
230 if (iter != end()){
231 ++iter;
232 atomIds.erase( key->getId() );
233 {
234 NOTIFY(AtomNrChanged);
235 atomIdPool.releaseId(key->getNr());
236 LocalToGlobalId.erase(key->getNr());
237 key->setNr(-1);
238 }
239 formula-=key->getType();
240 key->removeFromMolecule();
241 }
242 return iter;
243}
244
245pair<molecule::iterator,bool> molecule::insert ( atom * const key )
246{
247 OBSERVE;
248 NOTIFY(AtomInserted);
249 std::pair<iterator,bool> res = atomIds.insert(key->getId());
250 if (res.second) { // push atom if went well
251 NOTIFY(AtomNrChanged);
252 key->setNr(atomIdPool.getNextId());
253 InsertLocalToGlobalId(key);
254 setAtomName(key);
255 formula+=key->getType();
256 return res;
257 } else {
258 return pair<iterator,bool>(end(),res.second);
259 }
260}
261
262void molecule::setAtomName(atom *_atom) const
263{
264 std::stringstream sstr;
265 sstr << _atom->getType()->getSymbol() << _atom->getNr();
266 _atom->setName(sstr.str());
267}
268
269World::AtomComposite molecule::getAtomSet() const
270{
271 World::AtomComposite vector_of_atoms;
272 for (molecule::iterator iter = begin(); iter != end(); ++iter)
273 vector_of_atoms.push_back(*iter);
274 return vector_of_atoms;
275}
276
277/** Adds given atom \a *pointer from molecule list.
278 * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
279 * \param *pointer allocated and set atom
280 * \return true - succeeded, false - atom not found in list
281 */
282bool molecule::AddAtom(atom *pointer)
283{
284 if (pointer != NULL) {
285 // molecule::insert() is called by setMolecule()
286 pointer->setMolecule(this);
287 }
288 return true;
289};
290
291/** Adds a copy of the given atom \a *pointer from molecule list.
292 * Increases molecule::last_atom and gives last number to added atom.
293 * \param *pointer allocated and set atom
294 * \return pointer to the newly added atom
295 */
296atom * molecule::AddCopyAtom(atom *pointer)
297{
298 atom *retval = NULL;
299 if (pointer != NULL) {
300 atom *walker = pointer->clone();
301 AddAtom(walker);
302 retval=walker;
303 }
304 return retval;
305};
306
307/** Adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
308 * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
309 * a different scheme when adding \a *replacement atom for the given one.
310 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
311 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
312 * *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
313 * The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
314 * replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
315 * element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
316 * hydrogens forming this angle with *origin.
317 * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
318 * triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
319 * determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
320 * We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
321 * \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
322 * \f]
323 * vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
324 * the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
325 * The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
326 * the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
327 * \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
328 * \f]
329 * as the coordination of all three atoms in the coordinate system of these three vectors:
330 * \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
331 *
332 * \param *out output stream for debugging
333 * \param *Bond pointer to bond between \a *origin and \a *replacement
334 * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
335 * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
336 * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
337 * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
338 * \return number of atoms added, if < bond::BondDegree then something went wrong
339 * \todo double and triple bonds splitting (always use the tetraeder angle!)
340 */
341//bool molecule::AddHydrogenReplacementAtom(bond::ptr TopBond, atom *BottomOrigin, atom *TopOrigin, atom *TopReplacement, bool IsAngstroem)
342//{
343//// Info info(__func__);
344// bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
345// double bondlength; // bond length of the bond to be replaced/cut
346// double bondangle; // bond angle of the bond to be replaced/cut
347// double BondRescale; // rescale value for the hydrogen bond length
348// bond::ptr FirstBond;
349// bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
350// atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
351// double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
352// Vector Orthovector1, Orthovector2; // temporary vectors in coordination construction
353// Vector InBondvector; // vector in direction of *Bond
354// const RealSpaceMatrix &matrix = World::getInstance().getDomain().getM();
355// bond::ptr Binder;
356//
357// // create vector in direction of bond
358// InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
359// bondlength = InBondvector.Norm();
360//
361// // is greater than typical bond distance? Then we have to correct periodically
362// // the problem is not the H being out of the box, but InBondvector have the wrong direction
363// // due to TopReplacement or Origin being on the wrong side!
364// const BondGraph * const BG = World::getInstance().getBondGraph();
365// const range<double> MinMaxBondDistance(
366// BG->getMinMaxDistance(TopOrigin,TopReplacement));
367// if (!MinMaxBondDistance.isInRange(bondlength)) {
368//// LOG(4, "InBondvector is: " << InBondvector << ".");
369// Orthovector1.Zero();
370// for (int i=NDIM;i--;) {
371// l = TopReplacement->at(i) - TopOrigin->at(i);
372// if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
373// Orthovector1[i] = (l < 0) ? -1. : +1.;
374// } // (signs are correct, was tested!)
375// }
376// Orthovector1 *= matrix;
377// InBondvector -= Orthovector1; // subtract just the additional translation
378// bondlength = InBondvector.Norm();
379//// LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
380// } // periodic correction finished
381//
382// InBondvector.Normalize();
383// // get typical bond length and store as scale factor for later
384// ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
385// BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->getDegree()-1);
386// if (BondRescale == -1) {
387// ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
388// return false;
389// BondRescale = bondlength;
390// } else {
391// if (!IsAngstroem)
392// BondRescale /= (1.*AtomicLengthToAngstroem);
393// }
394//
395// // discern single, double and triple bonds
396// switch(TopBond->getDegree()) {
397// case 1:
398// FirstOtherAtom = World::getInstance().createAtom(); // new atom
399// FirstOtherAtom->setType(1); // element is Hydrogen
400// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
401// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
402// if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
403// FirstOtherAtom->father = TopReplacement;
404// BondRescale = bondlength;
405// } else {
406// FirstOtherAtom->father = NULL; // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
407// }
408// InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length
409// FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
410// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
411//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
412// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
413// Binder->Cyclic = false;
414// Binder->Type = GraphEdge::TreeEdge;
415// break;
416// case 2:
417// {
418// // determine two other bonds (warning if there are more than two other) plus valence sanity check
419// const BondList& ListOfBonds = TopOrigin->getListOfBonds();
420// for (BondList::const_iterator Runner = ListOfBonds.begin();
421// Runner != ListOfBonds.end();
422// ++Runner) {
423// if ((*Runner) != TopBond) {
424// if (FirstBond == NULL) {
425// FirstBond = (*Runner);
426// FirstOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
427// } else if (SecondBond == NULL) {
428// SecondBond = (*Runner);
429// SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin);
430// } else {
431// ELOG(2, "Detected more than four bonds for atom " << TopOrigin->getName());
432// }
433// }
434// }
435// }
436// if (SecondOtherAtom == NULL) { // then we have an atom with valence four, but only 3 bonds: one to replace and one which is TopBond (third is FirstBond)
437// SecondBond = TopBond;
438// SecondOtherAtom = TopReplacement;
439// }
440// if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
441//// LOG(3, "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane.");
442//
443// // determine the plane of these two with the *origin
444// try {
445// Orthovector1 = Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
446// }
447// catch(LinearDependenceException &excp){
448// LOG(0, boost::diagnostic_information(excp));
449// // TODO: figure out what to do with the Orthovector in this case
450// AllWentWell = false;
451// }
452// } else {
453// Orthovector1.GetOneNormalVector(InBondvector);
454// }
455// //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
456// // orthogonal vector and bond vector between origin and replacement form the new plane
457// Orthovector1.MakeNormalTo(InBondvector);
458// Orthovector1.Normalize();
459// //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
460//
461// // create the two Hydrogens ...
462// FirstOtherAtom = World::getInstance().createAtom();
463// SecondOtherAtom = World::getInstance().createAtom();
464// FirstOtherAtom->setType(1);
465// SecondOtherAtom->setType(1);
466// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
467// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
468// SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
469// SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
470// FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
471// SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
472// bondangle = TopOrigin->getType()->getHBondAngle(1);
473// if (bondangle == -1) {
474// ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->getDegree() << "!");
475// return false;
476// bondangle = 0;
477// }
478// bondangle *= M_PI/180./2.;
479//// LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
480//// LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
481// FirstOtherAtom->Zero();
482// SecondOtherAtom->Zero();
483// for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
484// FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
485// SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
486// }
487// FirstOtherAtom->Scale(BondRescale); // rescale by correct BondDistance
488// SecondOtherAtom->Scale(BondRescale);
489// //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
490// *FirstOtherAtom += TopOrigin->getPosition();
491// *SecondOtherAtom += TopOrigin->getPosition();
492// // ... and add to molecule
493// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
494// AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
495//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
496//// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
497// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
498// Binder->Cyclic = false;
499// Binder->Type = GraphEdge::TreeEdge;
500// Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
501// Binder->Cyclic = false;
502// Binder->Type = GraphEdge::TreeEdge;
503// break;
504// case 3:
505// // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
506// FirstOtherAtom = World::getInstance().createAtom();
507// SecondOtherAtom = World::getInstance().createAtom();
508// ThirdOtherAtom = World::getInstance().createAtom();
509// FirstOtherAtom->setType(1);
510// SecondOtherAtom->setType(1);
511// ThirdOtherAtom->setType(1);
512// FirstOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
513// FirstOtherAtom->setFixedIon(TopReplacement->getFixedIon());
514// SecondOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
515// SecondOtherAtom->setFixedIon(TopReplacement->getFixedIon());
516// ThirdOtherAtom->setAtomicVelocity(TopReplacement->getAtomicVelocity()); // copy velocity
517// ThirdOtherAtom->setFixedIon(TopReplacement->getFixedIon());
518// FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father
519// SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father
520// ThirdOtherAtom->father = NULL; // we are just an added hydrogen with no father
521//
522// // we need to vectors orthonormal the InBondvector
523// AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
524//// LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
525// try{
526// Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
527// }
528// catch(LinearDependenceException &excp) {
529// LOG(0, boost::diagnostic_information(excp));
530// AllWentWell = false;
531// }
532//// LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
533//
534// // create correct coordination for the three atoms
535// alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.; // retrieve triple bond angle from database
536// l = BondRescale; // desired bond length
537// b = 2.*l*sin(alpha); // base length of isosceles triangle
538// d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
539// f = b/sqrt(3.); // length for Orthvector1
540// g = b/2.; // length for Orthvector2
541//// LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
542//// LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
543// factors[0] = d;
544// factors[1] = f;
545// factors[2] = 0.;
546// FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
547// factors[1] = -0.5*f;
548// factors[2] = g;
549// SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
550// factors[2] = -g;
551// ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
552//
553// // rescale each to correct BondDistance
554//// FirstOtherAtom->x.Scale(&BondRescale);
555//// SecondOtherAtom->x.Scale(&BondRescale);
556//// ThirdOtherAtom->x.Scale(&BondRescale);
557//
558// // and relative to *origin atom
559// *FirstOtherAtom += TopOrigin->getPosition();
560// *SecondOtherAtom += TopOrigin->getPosition();
561// *ThirdOtherAtom += TopOrigin->getPosition();
562//
563// // ... and add to molecule
564// AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
565// AllWentWell = AllWentWell && AddAtom(SecondOtherAtom);
566// AllWentWell = AllWentWell && AddAtom(ThirdOtherAtom);
567//// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
568//// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
569//// LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
570// Binder = AddBond(BottomOrigin, FirstOtherAtom, 1);
571// Binder->Cyclic = false;
572// Binder->Type = GraphEdge::TreeEdge;
573// Binder = AddBond(BottomOrigin, SecondOtherAtom, 1);
574// Binder->Cyclic = false;
575// Binder->Type = GraphEdge::TreeEdge;
576// Binder = AddBond(BottomOrigin, ThirdOtherAtom, 1);
577// Binder->Cyclic = false;
578// Binder->Type = GraphEdge::TreeEdge;
579// break;
580// default:
581// ELOG(1, "BondDegree does not state single, double or triple bond!");
582// AllWentWell = false;
583// break;
584// }
585//
586// return AllWentWell;
587//};
588
589/** Creates a copy of this molecule.
590 * \param offset translation Vector for the new molecule relative to old one
591 * \return copy of molecule
592 */
593molecule *molecule::CopyMolecule(const Vector &offset) const
594{
595 molecule *copy = World::getInstance().createMolecule();
596
597 // copy all atoms
598 std::map< const atom *, atom *> FatherFinder;
599 for (iterator iter = begin(); iter != end(); ++iter) {
600 atom * const copy_atom = copy->AddCopyAtom(*iter);
601 copy_atom->setPosition(copy_atom->getPosition() + offset);
602 FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
603 }
604
605 // copy all bonds
606 for(const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
607 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
608 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
609 BondRunner != ListOfBonds.end();
610 ++BondRunner)
611 if ((*BondRunner)->leftatom == *AtomRunner) {
612 bond::ptr Binder = (*BondRunner);
613 // get the pendant atoms of current bond in the copy molecule
614 ASSERT(FatherFinder.count(Binder->leftatom),
615 "molecule::CopyMolecule() - No copy of original left atom "
616 +toString(Binder->leftatom)+" for bond copy found");
617 ASSERT(FatherFinder.count(Binder->rightatom),
618 "molecule::CopyMolecule() - No copy of original right atom "
619 +toString(Binder->rightatom)+" for bond copy found");
620 atom * const LeftAtom = FatherFinder[Binder->leftatom];
621 atom * const RightAtom = FatherFinder[Binder->rightatom];
622
623 bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
624 NewBond->Cyclic = Binder->Cyclic;
625 if (Binder->Cyclic)
626 copy->NoCyclicBonds++;
627 NewBond->Type = Binder->Type;
628 }
629 }
630 // correct fathers
631 //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
632
633 return copy;
634};
635
636
637/** Destroys all atoms inside this molecule.
638 */
639void molecule::removeAtomsinMolecule()
640{
641 // remove each atom from world
642 for(iterator AtomRunner = begin(); !empty(); AtomRunner = begin())
643 World::getInstance().destroyAtom(*AtomRunner);
644};
645
646
647/**
648 * Copies all atoms of a molecule which are within the defined parallelepiped.
649 *
650 * @param offest for the origin of the parallelepiped
651 * @param three vectors forming the matrix that defines the shape of the parallelpiped
652 */
653molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
654 molecule *copy = World::getInstance().createMolecule();
655
656 // copy all atoms
657 std::map< const atom *, atom *> FatherFinder;
658 for (iterator iter = begin(); iter != end(); ++iter) {
659 if (region.isInside((*iter)->getPosition())) {
660 atom * const copy_atom = copy->AddCopyAtom(*iter);
661 FatherFinder.insert( std::make_pair( *iter, copy_atom ) );
662 }
663 }
664
665 // copy all bonds
666 for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
667 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
668 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
669 BondRunner != ListOfBonds.end();
670 ++BondRunner)
671 if ((*BondRunner)->leftatom == *AtomRunner) {
672 bond::ptr Binder = (*BondRunner);
673 if ((FatherFinder.count(Binder->leftatom))
674 && (FatherFinder.count(Binder->rightatom))) {
675 // if copy present, then it must be from subregion
676 atom * const LeftAtom = FatherFinder[Binder->leftatom];
677 atom * const RightAtom = FatherFinder[Binder->rightatom];
678
679 bond::ptr const NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->getDegree());
680 NewBond->Cyclic = Binder->Cyclic;
681 if (Binder->Cyclic)
682 copy->NoCyclicBonds++;
683 NewBond->Type = Binder->Type;
684 }
685 }
686 }
687 // correct fathers
688 //for_each(begin(),end(),mem_fun(&atom::CorrectFather));
689
690 //TODO: copy->BuildInducedSubgraph(this);
691
692 return copy;
693}
694
695/** Adds a bond to a the molecule specified by two atoms, \a *first and \a *second.
696 * Also updates molecule::BondCount and molecule::NoNonBonds.
697 * \param *first first atom in bond
698 * \param *second atom in bond
699 * \return pointer to bond or NULL on failure
700 */
701bond::ptr molecule::AddBond(atom *atom1, atom *atom2, int degree)
702{
703 bond::ptr Binder;
704
705 // some checks to make sure we are able to create the bond
706 ASSERT(atom1,
707 "molecule::AddBond() - First atom "+toString(atom1)
708 +" is not a invalid pointer");
709 ASSERT(atom2,
710 "molecule::AddBond() - Second atom "+toString(atom2)
711 +" is not a invalid pointer");
712 ASSERT(isInMolecule(atom1),
713 "molecule::AddBond() - First atom "+toString(atom1)
714 +" is not part of molecule");
715 ASSERT(isInMolecule(atom2),
716 "molecule::AddBond() - Second atom "+toString(atom2)
717 +" is not part of molecule");
718
719 Binder.reset(new bond(atom1, atom2, degree));
720 atom1->RegisterBond(WorldTime::getTime(), Binder);
721 atom2->RegisterBond(WorldTime::getTime(), Binder);
722 if ((atom1->getType() != NULL)
723 && (atom1->getType()->getAtomicNumber() != 1)
724 && (atom2->getType() != NULL)
725 && (atom2->getType()->getAtomicNumber() != 1))
726 NoNonBonds++;
727
728 return Binder;
729};
730
731/** Set molecule::name from the basename without suffix in the given \a *filename.
732 * \param *filename filename
733 */
734void molecule::SetNameFromFilename(const char *filename)
735{
736 OBSERVE;
737 int length = 0;
738 const char *molname = strrchr(filename, '/');
739 if (molname != NULL)
740 molname += sizeof(char); // search for filename without dirs
741 else
742 molname = filename; // contains no slashes
743 const char *endname = strchr(molname, '.');
744 if ((endname == NULL) || (endname < molname))
745 length = strlen(molname);
746 else
747 length = strlen(molname) - strlen(endname);
748 cout << "Set name of molecule " << getId() << " to " << molname << endl;
749 strncpy(name, molname, length);
750 name[length]='\0';
751};
752
753/** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
754 * \param *dim vector class
755 */
756void molecule::SetBoxDimension(Vector *dim)
757{
758 RealSpaceMatrix domain;
759 for(int i =0; i<NDIM;++i)
760 domain.at(i,i) = dim->at(i);
761 World::getInstance().setDomain(domain);
762};
763
764/** Removes atom from molecule list, but does not delete it.
765 * \param *pointer atom to be removed
766 * \return true - succeeded, false - atom not found in list
767 */
768bool molecule::UnlinkAtom(atom *pointer)
769{
770 if (pointer == NULL)
771 return false;
772 pointer->removeFromMolecule();
773 return true;
774};
775
776/** Removes every atom from molecule list.
777 * \return true - succeeded, false - atom not found in list
778 */
779bool molecule::CleanupMolecule()
780{
781 for (molecule::iterator iter = begin(); !empty(); iter = begin())
782 (*iter)->removeFromMolecule();
783 return empty();
784};
785
786/** Finds an atom specified by its continuous number.
787 * \param Nr number of atom withim molecule
788 * \return pointer to atom or NULL
789 */
790atom * molecule::FindAtom(int Nr) const
791{
792 LocalToGlobalId_t::const_iterator iter = LocalToGlobalId.find(Nr);
793 if (iter != LocalToGlobalId.end()) {
794 //LOG(0, "Found Atom Nr. " << walker->getNr());
795 return iter->second;
796 } else {
797 ELOG(1, "Atom with Nr " << Nr << " not found in molecule " << getName() << "'s list.");
798 return NULL;
799 }
800}
801
802/** Checks whether the given atom is a member of this molecule.
803 *
804 * We make use here of molecule::atomIds to get a result on
805 *
806 * @param _atom atom to check
807 * @return true - is member, false - is not
808 */
809bool molecule::isInMolecule(const atom * const _atom)
810{
811 ASSERT(_atom->getMolecule() == this,
812 "molecule::isInMolecule() - atom is not designated to be in molecule '"
813 +toString(this->getName())+"'.");
814 molecule::const_iterator iter = atomIds.find(_atom->getId());
815 return (iter != atomIds.end());
816}
817
818/** Asks for atom number, and checks whether in list.
819 * \param *text question before entering
820 */
821atom * molecule::AskAtom(std::string text)
822{
823 int No;
824 atom *ion = NULL;
825 do {
826 //std::cout << "============Atom list==========================" << std::endl;
827 //mol->Output((ofstream *)&cout);
828 //std::cout << "===============================================" << std::endl;
829 std::cout << text;
830 cin >> No;
831 ion = this->FindAtom(No);
832 } while (ion == NULL);
833 return ion;
834};
835
836/** Checks if given coordinates are within cell volume.
837 * \param *x array of coordinates
838 * \return true - is within, false - out of cell
839 */
840bool molecule::CheckBounds(const Vector *x) const
841{
842 const RealSpaceMatrix &domain = World::getInstance().getDomain().getM();
843 bool result = true;
844 for (int i=0;i<NDIM;i++) {
845 result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
846 }
847 //return result;
848 return true; /// probably not gonna use the check no more
849};
850
851/** Prints molecule to *out.
852 * \param *out output stream
853 */
854bool molecule::Output(ostream * const output) const
855{
856 if (output == NULL) {
857 return false;
858 } else {
859 int AtomNo[MAX_ELEMENTS];
860 memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
861 enumeration<const element*> elementLookup = formula.enumerateElements();
862 *output << "#Ion_TypeNr._Nr.R[0] R[1] R[2] MoveType (0 MoveIon, 1 FixedIon)" << endl;
863 for_each(begin(),end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
864 return true;
865 }
866};
867
868/** Outputs contents of each atom::ListOfBonds.
869 * \param *out output stream
870 */
871void molecule::OutputListOfBonds() const
872{
873 std::stringstream output;
874 LOG(2, "From Contents of ListOfBonds, all atoms:");
875 for (molecule::const_iterator iter = begin();
876 iter != end();
877 ++iter) {
878 (*iter)->OutputBondOfAtom(output);
879 output << std::endl << "\t\t";
880 }
881 LOG(2, output.str());
882}
883
884/** Brings molecule::AtomCount and atom::*Name up-to-date.
885 * \param *out output stream for debugging
886 */
887size_t molecule::doCountNoNonHydrogen() const
888{
889 int temp = 0;
890 // go through atoms and look for new ones
891 for (molecule::const_iterator iter = begin(); iter != end(); ++iter)
892 if ((*iter)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it
893 ++temp;
894 return temp;
895};
896
897/** Counts the number of present bonds.
898 * \return number of bonds
899 */
900int molecule::doCountBonds() const
901{
902 unsigned int counter = 0;
903 for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
904 const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
905 for(BondList::const_iterator BondRunner = ListOfBonds.begin();
906 BondRunner != ListOfBonds.end();
907 ++BondRunner)
908 if ((*BondRunner)->leftatom == *AtomRunner)
909 counter++;
910 }
911 return counter;
912}
913
914
915/** Returns an index map for two father-son-molecules.
916 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
917 * \param *out output stream for debugging
918 * \param *OtherMolecule corresponding molecule with fathers
919 * \return allocated map of size molecule::AtomCount with map
920 * \todo make this with a good sort O(n), not O(n^2)
921 */
922int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
923{
924 LOG(3, "Begin of GetFatherAtomicMap.");
925 int *AtomicMap = new int[getAtomCount()];
926 for (int i=getAtomCount();i--;)
927 AtomicMap[i] = -1;
928 if (OtherMolecule == this) { // same molecule
929 for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
930 AtomicMap[i] = i;
931 LOG(4, "Map is trivial.");
932 } else {
933 std::stringstream output;
934 output << "Map is ";
935 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
936 if ((*iter)->father == NULL) {
937 AtomicMap[(*iter)->getNr()] = -2;
938 } else {
939 for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
940 //for (int i=0;i<AtomCount;i++) { // search atom
941 //for (int j=0;j<OtherMolecule->getAtomCount();j++) {
942 //LOG(4, "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << ".");
943 if ((*iter)->father == (*runner))
944 AtomicMap[(*iter)->getNr()] = (*runner)->getNr();
945 }
946 }
947 output << AtomicMap[(*iter)->getNr()] << "\t";
948 }
949 LOG(4, output.str());
950 }
951 LOG(3, "End of GetFatherAtomicMap.");
952 return AtomicMap;
953};
954
955
956void molecule::flipActiveFlag(){
957 ActiveFlag = !ActiveFlag;
958}
959
960Shape molecule::getBoundingShape(const double scale) const
961{
962 // create Sphere around every atom
963 if (empty())
964 return Nowhere();
965 const_iterator iter = begin();
966 const Vector center = (*iter)->getPosition();
967 const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
968 Shape BoundingShape = Sphere(center, vdWRadius*scale);
969 for(++iter; iter != end(); ++iter) {
970 const Vector center = (*iter)->getPosition();
971 const double vdWRadius = (*iter)->getElement().getVanDerWaalsRadius();
972 if (vdWRadius*scale != 0.)
973 BoundingShape = Sphere(center, vdWRadius*scale) || BoundingShape;
974 }
975 return BoundingShape;
976}
977
978Shape molecule::getBoundingSphere(const double boundary) const
979{
980 // get center and radius
981 Vector center;
982 double radius = 0.;
983 {
984 center.Zero();
985 for(const_iterator iter = begin(); iter != end(); ++iter)
986 center += (*iter)->getPosition();
987 center *= 1./(double)size();
988 for(const_iterator iter = begin(); iter != end(); ++iter) {
989 const Vector &position = (*iter)->getPosition();
990 const double temp_distance = position.DistanceSquared(center);
991 if (temp_distance > radius)
992 radius = temp_distance;
993 }
994 }
995 // convert radius to true value and add some small boundary
996 radius = sqrt(radius) + boundary + 1e+6*std::numeric_limits<double>::epsilon();
997 LOG(1, "INFO: The " << size() << " atoms of the molecule are contained in a sphere at "
998 << center << " with radius " << radius << ".");
999
1000 // TODO: When we do not use a Sphere here anymore, then FillRegularGridAction will
1001 // will not work as it expects a sphere due to possible random rotations.
1002 Shape BoundingShape(Sphere(center, radius));
1003 LOG(1, "INFO: Created sphere at " << BoundingShape.getCenter() << " and radius "
1004 << BoundingShape.getRadius() << ".");
1005 return BoundingShape;
1006}
1007
1008// construct idpool
1009CONSTRUCT_IDPOOL(atomId_t, continuousId)
1010
Note: See TracBrowser for help on using the repository browser.