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 |
|
---|
64 | SaturatedFragment::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 |
|
---|
92 | SaturatedFragment::~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 |
|
---|
110 | typedef std::vector<atom *> atoms_t;
|
---|
111 |
|
---|
112 | atoms_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 |
|
---|
128 | typedef std::map<atom *, BondList > CutBonds_t;
|
---|
129 |
|
---|
130 | CutBonds_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 |
|
---|
168 | typedef std::vector<atomId_t> atomids_t;
|
---|
169 |
|
---|
170 | atomids_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 |
|
---|
204 | void 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 |
|
---|
228 | void 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 |
|
---|
309 | const 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 |
|
---|
328 | bool 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 |
|
---|
346 | bool 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 |
|
---|
393 | atom * 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 |
|
---|
404 | bool 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
|
---|
621 | inline
|
---|
622 | #else
|
---|
623 | static
|
---|
624 | #endif
|
---|
625 | void 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 |
|
---|
632 | std::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 |
|
---|
672 | static 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 |
|
---|
711 | typedef std::map<atomId_t, atomId_t> idtable_t;
|
---|
712 |
|
---|
713 | static 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 |
|
---|
727 | static 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 |
|
---|
742 | static 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 |
|
---|
754 | FragmentationEdges::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 | }
|
---|