source: src/molecule.cpp@ c38826

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 c38826 was e4afb4, checked in by Frederik Heber <heber@…>, 14 years ago

Huge refactoring: Introduction of Traits to Actions.

This change is really big but the introduction of the Trait concept (at least
in its current light form) is so fundamental that lots of pieces had to be
changed in order to get everything working.

The main point why it was necessary to add these traits in the first place was
to comfortably allow for adding extension of Actions information-wise, i.e.
with stuff that is only important for the QtUI, such as icons, or tooltips, ...
This extra information should not be stored with Action itself, as it has
nothing to do with the workings of the Action. And neither should it get
stored with some blown-out-of-proportions MapOfActions class ...

The gist of the change is as follows:

  • OptionTrait contains the token, description, shortform and type of an option, such as ("position", "position in space, none, typeid(Vector)).
  • ActionTrait is the derived form for actions where additionally MenuPosition and MenuName are stored (and probably more to come for the GUI), also we have a set of OptionTrait instances, one for each option of the Action.
  • Action then contains this ActionTrait, specialized for each Action.
  • the preprocessor macros have been enhanced to gather all this information from the .def files.
  • MapOfActions is gone. Completely. Most of its use was to store this extra information and the ValueStorage part now is just in class ValueStorage.
  • ValueStorage is no more an interface to MapOfActions but as the name says a (type-safe) ValueStorage.

Listing the (remaining) changes in alphabetical order of the class:

  • Action
    • member value ::name dropped, ::getName() uses ActionTraits::getName()
    • new define NODEFAULT which is used in paramdefaults in .def files
    • all derived actions classes such as Process, Calculations, MakroAction,... have been adapated to use the ActionTrait concept as well.
  • ActionHistory
    • extraced RedoAction and UndoAction, shifted implementation into their own object files and they use .def files as well (i.e. streamlined with method used for other actions)
  • MenuDescription
    • contain information on Menus such as name, ...
    • new unit test checks for consistency
  • molecule
    • const member functions: Copy(), Output() and OutputBonds()
  • OptionRegistry
    • new registry class for options only
    • we want the same type throughout the code for each token, e.g. "position"
    • the registry containts checks for consistency
  • OptionTrait
    • default values are specified in paramdefaults, none are given by NODEFAULT
    • introduced default for translate-atoms, point-correlation, pair-correlation
  • Registry pattern
    • new unit test, but only sceleton code so far
  • ...Query, also ...Pipe
    • atoms, molecule and elements are now all const
    • also ValueStorage's signatures all have const therein
  • ValueStorage
    • set/queryCurrentValue from MapOfActions
    • at times VectorValue has been in .def files where Vector was in the signature. This is cleared. Such stuff is only present for e.g. BoxVector being queried as a Vector. But this is a feature and intended.
  • World
    • most of the (un)selection functions now work on const atoms and molecules
    • in one case we need a const_cast to remove this, but this is intentional, as the vector of selected atoms stores non-const pointers and this is ok.

There is only one test which had to be changed slightly because a specific
option token as "position" must now have the same type everywhere, e.g. always
Vector.

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