source: src/Fragmentation/Exporters/SaturatedFragment.cpp@ 9eb71b3

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes
Last change on this file since 9eb71b3 was 9eb71b3, checked in by Frederik Heber <frederik.heber@…>, 8 years ago

Commented out MemDebug include and Memory::ignore.

  • MemDebug clashes with various allocation operators that use a specific placement in memory. It is so far not possible to wrap new/delete fully. Hence, we stop this effort which so far has forced us to put ever more includes (with clashes) into MemDebug and thereby bloat compilation time.
  • MemDebug does not add that much usefulness which is not also provided by valgrind.
  • Property mode set to 100644
File size: 30.5 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2013 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * SaturatedFragment.cpp
26 *
27 * Created on: Mar 3, 2013
28 * Author: heber
29 */
30
31// include config.h
32#ifdef HAVE_CONFIG_H
33#include <config.h>
34#endif
35
36//#include "CodePatterns/MemDebug.hpp"
37
38#include "SaturatedFragment.hpp"
39
40#include <algorithm>
41#include <cmath>
42#include <iostream>
43
44#include <boost/bind.hpp>
45#include <boost/function.hpp>
46
47#include "CodePatterns/Assert.hpp"
48#include "CodePatterns/Log.hpp"
49
50#include "LinearAlgebra/Exceptions.hpp"
51#include "LinearAlgebra/Plane.hpp"
52#include "LinearAlgebra/RealSpaceMatrix.hpp"
53#include "LinearAlgebra/Vector.hpp"
54
55#include "Atom/atom.hpp"
56#include "Bond/bond.hpp"
57#include "config.hpp"
58#include "Descriptors/AtomIdDescriptor.hpp"
59#include "Fragmentation/Exporters/HydrogenPool.hpp"
60#include "Fragmentation/HydrogenSaturation_enum.hpp"
61#include "Graph/BondGraph.hpp"
62#include "World.hpp"
63
64SaturatedFragment::SaturatedFragment(
65 const KeySet &_set,
66 KeySetsInUse_t &_container,
67 HydrogenPool &_hydrogens,
68 const enum HydrogenTreatment _treatment,
69 const enum HydrogenSaturation _saturation,
70 const GlobalSaturationPositions_t &_globalsaturationpositions) :
71 container(_container),
72 set(_set),
73 hydrogens(_hydrogens),
74 FullMolecule(set),
75 treatment(_treatment),
76 saturation(_saturation)
77{
78 // add to in-use container
79 ASSERT( container.find(set) == container.end(),
80 "SaturatedFragment::SaturatedFragment() - the set "
81 +toString(set)+" is already marked as in use.");
82 container.insert(set);
83
84 // prepare saturation hydrogens, either using global information
85 // or if not given, local information (created in the function)
86 if (_globalsaturationpositions.empty())
87 saturate();
88 else
89 saturate(_globalsaturationpositions);
90}
91
92SaturatedFragment::~SaturatedFragment()
93{
94 // release all saturation hydrogens if present
95 for (KeySet::iterator iter = SaturationHydrogens.begin();
96 !SaturationHydrogens.empty();
97 iter = SaturationHydrogens.begin()) {
98 hydrogens.releaseHydrogen(*iter);
99 SaturationHydrogens.erase(iter);
100 }
101
102 // remove ourselves from in-use container
103 KeySetsInUse_t::iterator iter = container.find(set);
104 ASSERT( container.find(set) != container.end(),
105 "SaturatedFragment::SaturatedFragment() - the set "
106 +toString(set)+" is not marked as in use.");
107 container.erase(iter);
108}
109
110typedef std::vector<atom *> atoms_t;
111
112atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
113{
114 atoms_t atoms;
115 for (KeySet::const_iterator iter = _FullMolecule.begin();
116 iter != _FullMolecule.end();
117 ++iter) {
118 atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
119 ASSERT( Walker != NULL,
120 "gatherAllAtoms() - id "
121 +toString(*iter)+" is unknown to World.");
122 atoms.push_back(Walker);
123 }
124
125 return atoms;
126}
127
128typedef std::map<atom *, BondList > CutBonds_t;
129
130CutBonds_t gatherCutBonds(
131 const atoms_t &_atoms,
132 const KeySet &_set,
133 const enum HydrogenTreatment _treatment)
134{
135 // bool LonelyFlag = false;
136 CutBonds_t CutBonds;
137 for (atoms_t::const_iterator iter = _atoms.begin();
138 iter != _atoms.end();
139 ++iter) {
140 atom * const Walker = *iter;
141
142 // go through all bonds
143 const BondList& ListOfBonds = Walker->getListOfBonds();
144 for (BondList::const_iterator BondRunner = ListOfBonds.begin();
145 BondRunner != ListOfBonds.end();
146 ++BondRunner) {
147 atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
148 // if other atom is in key set or is a specially treated hydrogen
149 if (_set.find(OtherWalker->getId()) != _set.end()) {
150 LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
151 } else if ((_treatment == ExcludeHydrogen)
152 && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
153 LOG(4, "DEBUG: Walker " << *Walker << " is bound to specially treated hydrogen " <<
154 *OtherWalker << ".");
155 } else {
156 LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
157 << *OtherWalker << ", who is not in this fragment molecule.");
158 if (CutBonds.count(Walker) == 0)
159 CutBonds.insert( std::make_pair(Walker, BondList() ));
160 CutBonds[Walker].push_back(*BondRunner);
161 }
162 }
163 }
164
165 return CutBonds;
166}
167
168typedef std::vector<atomId_t> atomids_t;
169
170atomids_t gatherPresentExcludedHydrogens(
171 const atoms_t &_atoms,
172 const KeySet &_set,
173 const enum HydrogenTreatment _treatment)
174{
175 // bool LonelyFlag = false;
176 atomids_t ExcludedHydrogens;
177 for (atoms_t::const_iterator iter = _atoms.begin();
178 iter != _atoms.end();
179 ++iter) {
180 atom * const Walker = *iter;
181
182 // go through all bonds
183 const BondList& ListOfBonds = Walker->getListOfBonds();
184 for (BondList::const_iterator BondRunner = ListOfBonds.begin();
185 BondRunner != ListOfBonds.end();
186 ++BondRunner) {
187 atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
188 // if other atom is in key set or is a specially treated hydrogen
189 if (_set.find(OtherWalker->getId()) != _set.end()) {
190 LOG(6, "DEBUG: OtherWalker " << *OtherWalker << " is in set already.");
191 } else if ((_treatment == ExcludeHydrogen)
192 && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
193 LOG(5, "DEBUG: Adding excluded hydrogen OtherWalker " << *OtherWalker << ".");
194 ExcludedHydrogens.push_back(OtherWalker->getId());
195 } else {
196 LOG(6, "DEBUG: OtherWalker " << *Walker << " is not in this fragment molecule and no hydrogen.");
197 }
198 }
199 }
200
201 return ExcludedHydrogens;
202}
203
204void SaturatedFragment::saturate()
205{
206 // so far, we just have a set of keys. Hence, convert to atom refs
207 // and gather all atoms in a vector
208 std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
209
210 // go through each atom of the fragment and gather all cut bonds in list
211 CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
212
213 // add excluded hydrogens to FullMolecule if treated specially
214 if (treatment == ExcludeHydrogen) {
215 atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
216 FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
217 }
218
219 // go through all cut bonds and replace with a hydrogen
220 for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
221 atomiter != CutBonds.end(); ++atomiter) {
222 atom * const Walker = atomiter->first;
223 if (!saturateAtom(Walker, atomiter->second))
224 exit(1);
225 }
226}
227
228void SaturatedFragment::saturate(
229 const GlobalSaturationPositions_t &_globalsaturationpositions)
230{
231 // so far, we just have a set of keys. Hence, convert to atom refs
232 // and gather all atoms in a vector
233 std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
234
235 // go through each atom of the fragment and gather all cut bonds in list
236 CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
237
238 // add excluded hydrogens to FullMolecule if treated specially
239 if (treatment == ExcludeHydrogen) {
240 atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
241 FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
242 }
243
244 // go through all cut bonds and replace with a hydrogen
245 if (saturation == DoSaturate) {
246 for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
247 atomiter != CutBonds.end(); ++atomiter) {
248 atom * const Walker = atomiter->first;
249 ASSERT( set.find(Walker->getId()) != set.end(),
250 "SaturatedFragment::saturate() - key "
251 +toString(*Walker)+"not in set?");
252 LOG(4, "DEBUG: We are now saturating dangling bonds of " << *Walker);
253
254 // gather set of positions for this atom from global map
255 GlobalSaturationPositions_t::const_iterator mapiter =
256 _globalsaturationpositions.find(Walker->getId());
257 ASSERT( mapiter != _globalsaturationpositions.end(),
258 "SaturatedFragment::saturate() - no global information for "
259 +toString(*Walker));
260 const SaturationsPositionsPerNeighbor_t &saturationpositions =
261 mapiter->second;
262
263 // go through all cut bonds for this atom
264 for (BondList::const_iterator bonditer = atomiter->second.begin();
265 bonditer != atomiter->second.end(); ++bonditer) {
266 atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
267
268 // get positions from global map
269 SaturationsPositionsPerNeighbor_t::const_iterator positionsiter =
270 saturationpositions.find(OtherWalker->getId());
271 ASSERT(positionsiter != saturationpositions.end(),
272 "SaturatedFragment::saturate() - no information on bond neighbor "
273 +toString(*OtherWalker)+" to atom "+toString(*Walker));
274 ASSERT(!positionsiter->second.empty(),
275 "SaturatedFragment::saturate() - no positions for saturating bond to"
276 +toString(*OtherWalker)+" to atom "+toString(*Walker));
277
278// // get typical bond distance from elements database
279// double BondDistance = Walker->getType()->getHBondDistance(positionsiter->second.size()-1);
280// if (BondDistance < 0.) {
281// ELOG(2, "saturateAtoms() - no typical hydrogen bond distance of degree "
282// +toString(positionsiter->second.size())+" for element "
283// +toString(Walker->getType()->getName()));
284// // try bond degree 1 distance
285// BondDistance = Walker->getType()->getHBondDistance(1-1);
286// if (BondDistance < 0.) {
287// ELOG(1, "saturateAtoms() - no typical hydrogen bond distance for element "
288// +toString(Walker->getType()->getName()));
289// BondDistance = 1.;
290// }
291// }
292
293 // place hydrogen at each point
294 LOG(4, "DEBUG: Places to saturate for atom " << *OtherWalker
295 << " are " << positionsiter->second);
296 atom * const father = OtherWalker;
297 for (SaturationsPositions_t::const_iterator positer = positionsiter->second.begin();
298 positer != positionsiter->second.end(); ++positer) {
299 const atom& hydrogen =
300 setHydrogenReplacement(Walker, *positer, 1., father);
301 FullMolecule.insert(hydrogen.getId());
302 }
303 }
304 }
305 } else
306 LOG(3, "INFO: We are not saturating cut bonds.");
307}
308
309const atom& SaturatedFragment::setHydrogenReplacement(
310 const atom * const _OwnerAtom,
311 const Vector &_position,
312 const double _distance,
313 atom * const _father)
314{
315 atom * const _atom = hydrogens.leaseHydrogen(); // new atom
316 _atom->setPosition( _OwnerAtom->getPosition() + _distance * _position );
317 // always set as fixed ion (not moving during molecular dynamics simulation)
318 _atom->setFixedIon(true);
319 // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
320 _atom->setFather(_father);
321 SaturationHydrogens.insert(_atom->getId());
322 // note down the id of the replaced atom
323 replaced_atoms.insert( std::make_pair(_father->getId(), _atom->getId()) );
324
325 return *_atom;
326}
327
328bool SaturatedFragment::saturateAtom(
329 atom * const _atom,
330 const BondList &_cutbonds)
331{
332 // go through each bond and replace
333 for (BondList::const_iterator bonditer = _cutbonds.begin();
334 bonditer != _cutbonds.end(); ++bonditer) {
335 atom * const OtherWalker = (*bonditer)->GetOtherAtom(_atom);
336 if (!AddHydrogenReplacementAtom(
337 (*bonditer),
338 _atom,
339 OtherWalker,
340 World::getInstance().getConfig()->IsAngstroem == 1))
341 return false;
342 }
343 return true;
344}
345
346bool SaturatedFragment::OutputConfig(
347 std::ostream &out,
348 const ParserTypes _type) const
349{
350 // gather all atoms in a vector
351 std::vector<const atom *> atoms;
352 for (KeySet::const_iterator iter = FullMolecule.begin();
353 iter != FullMolecule.end();
354 ++iter) {
355 const atom * const Walker = const_cast<const World &>(World::getInstance()).
356 getAtom(AtomById(*iter));
357 ASSERT( Walker != NULL,
358 "SaturatedFragment::OutputConfig() - id "
359 +toString(*iter)+" is unknown to World.");
360 atoms.push_back(Walker);
361 }
362
363 // TODO: ScanForPeriodicCorrection() is missing so far!
364 // note however that this is not straight-forward for the following reasons:
365 // - we do not copy all atoms anymore, hence we are forced to shift the real
366 // atoms hither and back again
367 // - we use a long-range potential that supports periodic boundary conditions.
368 // Hence, there we would like the original configuration (split through the
369 // the periodic boundaries). Otherwise, we would have to shift (and probably
370 // interpolate) the potential with OBCs applying.
371
372 // list atoms in fragment for debugging
373 {
374 std::stringstream output;
375 output << "INFO: Contained atoms: ";
376 for (std::vector<const atom *>::const_iterator iter = atoms.begin();
377 iter != atoms.end(); ++iter) {
378 output << (*iter)->getName() << " ";
379 }
380 LOG(3, output.str());
381 }
382
383 // store to stream via FragmentParser
384 const bool intermediateResult =
385 FormatParserStorage::getInstance().save(
386 out,
387 FormatParserStorage::getInstance().getSuffixFromType(_type),
388 atoms);
389
390 return intermediateResult;
391}
392
393atom * const SaturatedFragment::getHydrogenReplacement(atom * const replacement)
394{
395 atom * const _atom = hydrogens.leaseHydrogen(); // new atom
396 _atom->setAtomicVelocity(replacement->getAtomicVelocity()); // copy velocity
397 _atom->setFixedIon(replacement->getFixedIon());
398 // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
399 _atom->setFather(replacement);
400 SaturationHydrogens.insert(_atom->getId());
401 return _atom;
402}
403
404bool SaturatedFragment::AddHydrogenReplacementAtom(
405 bond::ptr TopBond,
406 atom *Origin,
407 atom *Replacement,
408 bool IsAngstroem)
409{
410// Info info(__func__);
411 bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
412 double bondlength; // bond length of the bond to be replaced/cut
413 double bondangle; // bond angle of the bond to be replaced/cut
414 double BondRescale; // rescale value for the hydrogen bond length
415 bond::ptr FirstBond;
416 bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
417 atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
418 double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
419 Vector Orthovector1, Orthovector2; // temporary vectors in coordination construction
420 Vector InBondvector; // vector in direction of *Bond
421 const RealSpaceMatrix &matrix = World::getInstance().getDomain().getM();
422 bond::ptr Binder;
423
424 // create vector in direction of bond
425 InBondvector = Replacement->getPosition() - Origin->getPosition();
426 bondlength = InBondvector.Norm();
427
428 // is greater than typical bond distance? Then we have to correct periodically
429 // the problem is not the H being out of the box, but InBondvector have the wrong direction
430 // due to Replacement or Origin being on the wrong side!
431 const BondGraph * const BG = World::getInstance().getBondGraph();
432 const range<double> MinMaxBondDistance(
433 BG->getMinMaxDistance(Origin,Replacement));
434 if (!MinMaxBondDistance.isInRange(bondlength)) {
435// LOG(4, "InBondvector is: " << InBondvector << ".");
436 Orthovector1.Zero();
437 for (int i=NDIM;i--;) {
438 l = Replacement->at(i) - Origin->at(i);
439 if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
440 Orthovector1[i] = (l < 0) ? -1. : +1.;
441 } // (signs are correct, was tested!)
442 }
443 Orthovector1 *= matrix;
444 InBondvector -= Orthovector1; // subtract just the additional translation
445 bondlength = InBondvector.Norm();
446// LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
447 } // periodic correction finished
448
449 InBondvector.Normalize();
450 // get typical bond length and store as scale factor for later
451 ASSERT(Origin->getType() != NULL,
452 "SaturatedFragment::AddHydrogenReplacementAtom() - element of Origin is not given.");
453 BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree()-1);
454 if (BondRescale == -1) {
455 ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
456 BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
457 if (BondRescale == -1) {
458 ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
459 return false;
460 BondRescale = bondlength;
461 }
462 } else {
463 if (!IsAngstroem)
464 BondRescale /= (1.*AtomicLengthToAngstroem);
465 }
466
467 // discern single, double and triple bonds
468 switch(TopBond->getDegree()) {
469 case 1:
470 // check whether replacement has been an excluded hydrogen
471 if (Replacement->getType()->getAtomicNumber() == HydrogenPool::HYDROGEN) { // neither rescale nor replace if it's already hydrogen
472 FirstOtherAtom = Replacement;
473 BondRescale = bondlength;
474 } else {
475 FirstOtherAtom = getHydrogenReplacement(Replacement);
476 InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length
477 FirstOtherAtom->setPosition(Origin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
478 }
479 FullMolecule.insert(FirstOtherAtom->getId());
480// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
481 break;
482 case 2:
483 {
484 // determine two other bonds (warning if there are more than two other) plus valence sanity check
485 const BondList& ListOfBonds = Origin->getListOfBonds();
486 for (BondList::const_iterator Runner = ListOfBonds.begin();
487 Runner != ListOfBonds.end();
488 ++Runner) {
489 if ((*Runner) != TopBond) {
490 if (FirstBond == NULL) {
491 FirstBond = (*Runner);
492 FirstOtherAtom = (*Runner)->GetOtherAtom(Origin);
493 } else if (SecondBond == NULL) {
494 SecondBond = (*Runner);
495 SecondOtherAtom = (*Runner)->GetOtherAtom(Origin);
496 } else {
497 ELOG(2, "Detected more than four bonds for atom " << Origin->getName());
498 }
499 }
500 }
501 }
502 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)
503 SecondBond = TopBond;
504 SecondOtherAtom = Replacement;
505 }
506 if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
507// LOG(3, "Regarding the double bond (" << Origin->Name << "<->" << Replacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << Origin->Name << " to determine orthogonal plane.");
508
509 // determine the plane of these two with the *origin
510 try {
511 Orthovector1 = Plane(Origin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
512 }
513 catch(LinearDependenceException &excp){
514 LOG(0, boost::diagnostic_information(excp));
515 // TODO: figure out what to do with the Orthovector in this case
516 AllWentWell = false;
517 }
518 } else {
519 Orthovector1.GetOneNormalVector(InBondvector);
520 }
521 //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
522 // orthogonal vector and bond vector between origin and replacement form the new plane
523 Orthovector1.MakeNormalTo(InBondvector);
524 Orthovector1.Normalize();
525 //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
526
527 // create the two Hydrogens ...
528 FirstOtherAtom = getHydrogenReplacement(Replacement);
529 SecondOtherAtom = getHydrogenReplacement(Replacement);
530 FullMolecule.insert(FirstOtherAtom->getId());
531 FullMolecule.insert(SecondOtherAtom->getId());
532 bondangle = Origin->getType()->getHBondAngle(1);
533 if (bondangle == -1) {
534 ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
535 return false;
536 bondangle = 0;
537 }
538 bondangle *= M_PI/180./2.;
539// LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
540// LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
541 FirstOtherAtom->Zero();
542 SecondOtherAtom->Zero();
543 for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
544 FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
545 SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
546 }
547 FirstOtherAtom->Scale(BondRescale); // rescale by correct BondDistance
548 SecondOtherAtom->Scale(BondRescale);
549 //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
550 *FirstOtherAtom += Origin->getPosition();
551 *SecondOtherAtom += Origin->getPosition();
552 // ... and add to molecule
553// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
554// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
555 break;
556 case 3:
557 // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
558 FirstOtherAtom = getHydrogenReplacement(Replacement);
559 SecondOtherAtom = getHydrogenReplacement(Replacement);
560 ThirdOtherAtom = getHydrogenReplacement(Replacement);
561 FullMolecule.insert(FirstOtherAtom->getId());
562 FullMolecule.insert(SecondOtherAtom->getId());
563 FullMolecule.insert(ThirdOtherAtom->getId());
564
565 // we need to vectors orthonormal the InBondvector
566 AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
567// LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
568 try{
569 Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
570 }
571 catch(LinearDependenceException &excp) {
572 LOG(0, boost::diagnostic_information(excp));
573 AllWentWell = false;
574 }
575// LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
576
577 // create correct coordination for the three atoms
578 alpha = (Origin->getType()->getHBondAngle(2))/180.*M_PI/2.; // retrieve triple bond angle from database
579 l = BondRescale; // desired bond length
580 b = 2.*l*sin(alpha); // base length of isosceles triangle
581 d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
582 f = b/sqrt(3.); // length for Orthvector1
583 g = b/2.; // length for Orthvector2
584// LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
585// LOG(3, "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g));
586 factors[0] = d;
587 factors[1] = f;
588 factors[2] = 0.;
589 FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
590 factors[1] = -0.5*f;
591 factors[2] = g;
592 SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
593 factors[2] = -g;
594 ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
595
596 // rescale each to correct BondDistance
597// FirstOtherAtom->x.Scale(&BondRescale);
598// SecondOtherAtom->x.Scale(&BondRescale);
599// ThirdOtherAtom->x.Scale(&BondRescale);
600
601 // and relative to *origin atom
602 *FirstOtherAtom += Origin->getPosition();
603 *SecondOtherAtom += Origin->getPosition();
604 *ThirdOtherAtom += Origin->getPosition();
605
606 // ... and add to molecule
607// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
608// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
609// LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
610 break;
611 default:
612 ELOG(1, "BondDegree does not state single, double or triple bond!");
613 AllWentWell = false;
614 break;
615 }
616
617 return AllWentWell;
618};
619
620#ifdef HAVE_INLINE
621inline
622#else
623static
624#endif
625void updateVector(Vector &_first, const Vector &_second,
626 const boost::function<const double& (const double &, const double &)> &_comparator)
627{
628 for (size_t i=0;i<NDIM;++i)
629 _first[i] = _comparator(_first[i], _second[i]);
630}
631
632std::pair<Vector, Vector> SaturatedFragment::getRoughBoundingBox() const
633{
634 typedef boost::function<const double& (const double &, const double &)> comparator_t;
635 static const comparator_t minimizer = boost::bind(&std::min<double>, _1, _2);
636 static const comparator_t maximizer = boost::bind(&std::max<double>, _1, _2);
637 // initialize return values
638 Vector minimum;
639 Vector maximum;
640 for (size_t i=0;i<NDIM;++i) {
641 minimum[i] = std::numeric_limits<double>::max();
642 maximum[i] = -std::numeric_limits<double>::max();
643 }
644
645 // go through all "core" atoms of the fragment
646 const std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
647 for (std::vector<atom *>::const_iterator iter = atoms.begin();
648 iter != atoms.end(); ++iter) {
649 const Vector &position = (*iter)->getPosition();
650 updateVector(minimum, position, minimizer);
651 updateVector(maximum, position, maximizer);
652 }
653
654 // go through each atom of the fragment and gather all cut bonds in list
655 const CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
656 for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
657 atomiter != CutBonds.end(); ++atomiter) {
658 const atom * const walker = atomiter->first;
659 const BondList &cutbonds = atomiter->second;
660 for (BondList::const_iterator bonditer = cutbonds.begin();
661 bonditer != cutbonds.end(); ++bonditer) {
662 const atom * const OtherWalker = (*bonditer)->GetOtherAtom(walker);
663 const Vector &position = OtherWalker->getPosition();
664 updateVector(minimum, position, minimizer);
665 updateVector(maximum, position, maximizer);
666 }
667 }
668
669 return std::make_pair(minimum, maximum);
670}
671
672static FragmentationEdges::edges_t createEdgesFromAtoms(
673 const std::vector<atom *> &_atoms,
674 const SaturatedFragment::replaced_atoms_t &replaced_atoms,
675 const KeySet &_FullMolecule)
676{
677 FragmentationEdges::edges_t edges;
678 for (std::vector<atom *>::const_iterator iter = _atoms.begin();
679 iter != _atoms.end(); ++iter) {
680 const atom * const walker = *iter;
681 const atomId_t &walkerid = walker->getId();
682 const BondList &ListOfBonds = walker->getListOfBonds();
683 for (BondList::const_iterator bonditer = ListOfBonds.begin();
684 bonditer != ListOfBonds.end(); ++bonditer) {
685 const atom * const OtherWalker = (*bonditer)->GetOtherAtom(walker);
686 const atomId_t &otherwalkerid = OtherWalker->getId();
687 if (_FullMolecule.count(otherwalkerid)) {
688 if (walkerid < otherwalkerid) {
689 // check if it's in fullkeysets (also contains excluded and saturation hydrogens)
690 if (_FullMolecule.count(otherwalkerid))
691 edges.push_back( FragmentationEdges::edge_t(walkerid, otherwalkerid) );
692 }
693 } else {
694 ASSERT( replaced_atoms.count(otherwalkerid),
695 "createEdgesFromAtoms() - atom #"+toString(otherwalkerid)
696 +" to a cut bond is not in replaced_atoms.");
697 // add bond to every saturation hydrogen instead
698 const std::pair<
699 SaturatedFragment::replaced_atoms_t::const_iterator,
700 SaturatedFragment::replaced_atoms_t::const_iterator> range = replaced_atoms.equal_range(otherwalkerid);
701 for (SaturatedFragment::replaced_atoms_t::const_iterator replaceiter = range.first;
702 replaceiter != range.second; ++replaceiter) {
703 edges.push_back( FragmentationEdges::edge_t(walkerid, replaceiter->second) );
704 }
705 }
706 }
707 }
708 return edges;
709}
710
711typedef std::map<atomId_t, atomId_t> idtable_t;
712
713static idtable_t createIdTable(
714 const std::vector<atom *> &_atoms)
715{
716 idtable_t idtable;
717 atomId_t newid = 0;
718 for (std::vector<atom *>::const_iterator iter = _atoms.begin();
719 iter != _atoms.end(); ++iter) {
720 const atom * const walker = *iter;
721 const atomId_t &walkerid = walker->getId();
722 idtable.insert( std::make_pair(walkerid, newid++) );
723 }
724 return idtable;
725}
726
727static atomId_t getTranslatedId(
728 const idtable_t &_idtable,
729 const atomId_t &_id
730 )
731{
732 idtable_t::const_iterator iter = _idtable.find(_id);
733 if (iter != _idtable.end())
734 return iter->second;
735 else {
736 ASSERT( 0,
737 "getTranslatedId() - idtable does not contain id in FullMolecule?");
738 return (atomId_t)-1;
739 }
740}
741
742static void rewriteEdgeIndices(
743 FragmentationEdges::edges_t &_edges,
744 const idtable_t &_idtable)
745{
746 for (FragmentationEdges::edges_t::iterator edgeiter = _edges.begin();
747 edgeiter != _edges.end(); ++edgeiter) {
748 FragmentationEdges::edge_t &edge = *edgeiter;
749 edge.first = getTranslatedId(_idtable, edge.first);
750 edge.second = getTranslatedId(_idtable, edge.second);
751 }
752}
753
754FragmentationEdges::edges_t SaturatedFragment::getEdges() const
755{
756 // gather all edges
757 const std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
758 FragmentationEdges::edges_t edges = createEdgesFromAtoms(atoms, replaced_atoms, FullMolecule);
759
760 // translate each edge
761 const std::map<atomId_t, atomId_t> idtable = createIdTable(atoms);
762
763 // rewrite indices of edges in correct order
764 rewriteEdgeIndices(edges, idtable);
765
766 // for debugging output edges
767 LOG(2, "DEBUG: FullMolecule is : " << FullMolecule << " with edges " << edges);
768
769 return edges;
770}
Note: See TracBrowser for help on using the repository browser.