source: src/molecule.cpp@ 4ef9b7

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 4ef9b7 was aafd77, checked in by Tillmann Crueger <crueger@…>, 14 years ago

Removed all inclusions of GSL-Headers from molecule.hpp

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