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 <cmath>
|
---|
41 | #include <iostream>
|
---|
42 |
|
---|
43 | #include "CodePatterns/Assert.hpp"
|
---|
44 | #include "CodePatterns/Log.hpp"
|
---|
45 |
|
---|
46 | #include "LinearAlgebra/Exceptions.hpp"
|
---|
47 | #include "LinearAlgebra/Plane.hpp"
|
---|
48 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
49 | #include "LinearAlgebra/Vector.hpp"
|
---|
50 |
|
---|
51 | #include "Atom/atom.hpp"
|
---|
52 | #include "Bond/bond.hpp"
|
---|
53 | #include "config.hpp"
|
---|
54 | #include "Descriptors/AtomIdDescriptor.hpp"
|
---|
55 | #include "Fragmentation/Exporters/HydrogenPool.hpp"
|
---|
56 | #include "Fragmentation/HydrogenSaturation_enum.hpp"
|
---|
57 | #include "Graph/BondGraph.hpp"
|
---|
58 | #include "World.hpp"
|
---|
59 |
|
---|
60 | SaturatedFragment::SaturatedFragment(
|
---|
61 | const KeySet &_set,
|
---|
62 | KeySetsInUse_t &_container,
|
---|
63 | HydrogenPool &_hydrogens,
|
---|
64 | const enum HydrogenTreatment _treatment,
|
---|
65 | const enum HydrogenSaturation _saturation,
|
---|
66 | const GlobalSaturationPositions_t &_globalsaturationpositions) :
|
---|
67 | container(_container),
|
---|
68 | set(_set),
|
---|
69 | hydrogens(_hydrogens),
|
---|
70 | FullMolecule(set),
|
---|
71 | treatment(_treatment),
|
---|
72 | saturation(_saturation)
|
---|
73 | {
|
---|
74 | // add to in-use container
|
---|
75 | ASSERT( container.find(set) == container.end(),
|
---|
76 | "SaturatedFragment::SaturatedFragment() - the set "
|
---|
77 | +toString(set)+" is already marked as in use.");
|
---|
78 | container.insert(set);
|
---|
79 |
|
---|
80 | // prepare saturation hydrogens, either using global information
|
---|
81 | // or if not given, local information (created in the function)
|
---|
82 | if (_globalsaturationpositions.empty())
|
---|
83 | saturate();
|
---|
84 | else
|
---|
85 | saturate(_globalsaturationpositions);
|
---|
86 | }
|
---|
87 |
|
---|
88 | SaturatedFragment::~SaturatedFragment()
|
---|
89 | {
|
---|
90 | // release all saturation hydrogens if present
|
---|
91 | for (KeySet::iterator iter = SaturationHydrogens.begin();
|
---|
92 | !SaturationHydrogens.empty();
|
---|
93 | iter = SaturationHydrogens.begin()) {
|
---|
94 | hydrogens.releaseHydrogen(*iter);
|
---|
95 | SaturationHydrogens.erase(iter);
|
---|
96 | }
|
---|
97 |
|
---|
98 | // remove ourselves from in-use container
|
---|
99 | KeySetsInUse_t::iterator iter = container.find(set);
|
---|
100 | ASSERT( container.find(set) != container.end(),
|
---|
101 | "SaturatedFragment::SaturatedFragment() - the set "
|
---|
102 | +toString(set)+" is not marked as in use.");
|
---|
103 | container.erase(iter);
|
---|
104 | }
|
---|
105 |
|
---|
106 | typedef std::vector<atom *> atoms_t;
|
---|
107 |
|
---|
108 | atoms_t gatherAllAtoms(const KeySet &_FullMolecule)
|
---|
109 | {
|
---|
110 | atoms_t atoms;
|
---|
111 | for (KeySet::const_iterator iter = _FullMolecule.begin();
|
---|
112 | iter != _FullMolecule.end();
|
---|
113 | ++iter) {
|
---|
114 | atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
|
---|
115 | ASSERT( Walker != NULL,
|
---|
116 | "gatherAllAtoms() - id "
|
---|
117 | +toString(*iter)+" is unknown to World.");
|
---|
118 | atoms.push_back(Walker);
|
---|
119 | }
|
---|
120 |
|
---|
121 | return atoms;
|
---|
122 | }
|
---|
123 |
|
---|
124 | typedef std::map<atom *, BondList > CutBonds_t;
|
---|
125 |
|
---|
126 | CutBonds_t gatherCutBonds(
|
---|
127 | const atoms_t &_atoms,
|
---|
128 | const KeySet &_set,
|
---|
129 | const enum HydrogenTreatment _treatment)
|
---|
130 | {
|
---|
131 | // bool LonelyFlag = false;
|
---|
132 | CutBonds_t CutBonds;
|
---|
133 | for (atoms_t::const_iterator iter = _atoms.begin();
|
---|
134 | iter != _atoms.end();
|
---|
135 | ++iter) {
|
---|
136 | atom * const Walker = *iter;
|
---|
137 |
|
---|
138 | // go through all bonds
|
---|
139 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
140 | for (BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
141 | BondRunner != ListOfBonds.end();
|
---|
142 | ++BondRunner) {
|
---|
143 | atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
|
---|
144 | // if other atom is in key set or is a specially treated hydrogen
|
---|
145 | if (_set.find(OtherWalker->getId()) != _set.end()) {
|
---|
146 | LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
|
---|
147 | } else if ((_treatment == ExcludeHydrogen)
|
---|
148 | && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
|
---|
149 | LOG(4, "DEBUG: Walker " << *Walker << " is bound to specially treated hydrogen " <<
|
---|
150 | *OtherWalker << ".");
|
---|
151 | } else {
|
---|
152 | LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
|
---|
153 | << *OtherWalker << ", who is not in this fragment molecule.");
|
---|
154 | if (CutBonds.count(Walker) == 0)
|
---|
155 | CutBonds.insert( std::make_pair(Walker, BondList() ));
|
---|
156 | CutBonds[Walker].push_back(*BondRunner);
|
---|
157 | }
|
---|
158 | }
|
---|
159 | }
|
---|
160 |
|
---|
161 | return CutBonds;
|
---|
162 | }
|
---|
163 |
|
---|
164 | typedef std::vector<atomId_t> atomids_t;
|
---|
165 |
|
---|
166 | atomids_t gatherPresentExcludedHydrogens(
|
---|
167 | const atoms_t &_atoms,
|
---|
168 | const KeySet &_set,
|
---|
169 | const enum HydrogenTreatment _treatment)
|
---|
170 | {
|
---|
171 | // bool LonelyFlag = false;
|
---|
172 | atomids_t ExcludedHydrogens;
|
---|
173 | for (atoms_t::const_iterator iter = _atoms.begin();
|
---|
174 | iter != _atoms.end();
|
---|
175 | ++iter) {
|
---|
176 | atom * const Walker = *iter;
|
---|
177 |
|
---|
178 | // go through all bonds
|
---|
179 | const BondList& ListOfBonds = Walker->getListOfBonds();
|
---|
180 | for (BondList::const_iterator BondRunner = ListOfBonds.begin();
|
---|
181 | BondRunner != ListOfBonds.end();
|
---|
182 | ++BondRunner) {
|
---|
183 | atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
|
---|
184 | // if other atom is in key set or is a specially treated hydrogen
|
---|
185 | if (_set.find(OtherWalker->getId()) != _set.end()) {
|
---|
186 | LOG(6, "DEBUG: OtherWalker " << *OtherWalker << " is in set already.");
|
---|
187 | } else if ((_treatment == ExcludeHydrogen)
|
---|
188 | && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
|
---|
189 | LOG(5, "DEBUG: Adding excluded hydrogen OtherWalker " << *OtherWalker << ".");
|
---|
190 | ExcludedHydrogens.push_back(OtherWalker->getId());
|
---|
191 | } else {
|
---|
192 | LOG(6, "DEBUG: OtherWalker " << *Walker << " is not in this fragment molecule and no hydrogen.");
|
---|
193 | }
|
---|
194 | }
|
---|
195 | }
|
---|
196 |
|
---|
197 | return ExcludedHydrogens;
|
---|
198 | }
|
---|
199 |
|
---|
200 | void SaturatedFragment::saturate()
|
---|
201 | {
|
---|
202 | // so far, we just have a set of keys. Hence, convert to atom refs
|
---|
203 | // and gather all atoms in a vector
|
---|
204 | std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
|
---|
205 |
|
---|
206 | // go through each atom of the fragment and gather all cut bonds in list
|
---|
207 | CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
|
---|
208 |
|
---|
209 | // add excluded hydrogens to FullMolecule if treated specially
|
---|
210 | if (treatment == ExcludeHydrogen) {
|
---|
211 | atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
|
---|
212 | FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
|
---|
213 | }
|
---|
214 |
|
---|
215 | // go through all cut bonds and replace with a hydrogen
|
---|
216 | for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
|
---|
217 | atomiter != CutBonds.end(); ++atomiter) {
|
---|
218 | atom * const Walker = atomiter->first;
|
---|
219 | if (!saturateAtom(Walker, atomiter->second))
|
---|
220 | exit(1);
|
---|
221 | }
|
---|
222 | }
|
---|
223 |
|
---|
224 | void SaturatedFragment::saturate(
|
---|
225 | const GlobalSaturationPositions_t &_globalsaturationpositions)
|
---|
226 | {
|
---|
227 | // so far, we just have a set of keys. Hence, convert to atom refs
|
---|
228 | // and gather all atoms in a vector
|
---|
229 | std::vector<atom *> atoms = gatherAllAtoms(FullMolecule);
|
---|
230 |
|
---|
231 | // go through each atom of the fragment and gather all cut bonds in list
|
---|
232 | CutBonds_t CutBonds = gatherCutBonds(atoms, set, treatment);
|
---|
233 |
|
---|
234 | // add excluded hydrogens to FullMolecule if treated specially
|
---|
235 | if (treatment == ExcludeHydrogen) {
|
---|
236 | atomids_t ExcludedHydrogens = gatherPresentExcludedHydrogens(atoms, set, treatment);
|
---|
237 | FullMolecule.insert(ExcludedHydrogens.begin(), ExcludedHydrogens.end());
|
---|
238 | }
|
---|
239 |
|
---|
240 | // go through all cut bonds and replace with a hydrogen
|
---|
241 | if (saturation == DoSaturate) {
|
---|
242 | for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
|
---|
243 | atomiter != CutBonds.end(); ++atomiter) {
|
---|
244 | atom * const Walker = atomiter->first;
|
---|
245 | LOG(4, "DEBUG: We are now saturating dangling bonds of " << *Walker);
|
---|
246 |
|
---|
247 | // gather set of positions for this atom from global map
|
---|
248 | GlobalSaturationPositions_t::const_iterator mapiter =
|
---|
249 | _globalsaturationpositions.find(Walker->getId());
|
---|
250 | ASSERT( mapiter != _globalsaturationpositions.end(),
|
---|
251 | "SaturatedFragment::saturate() - no global information for "
|
---|
252 | +toString(*Walker));
|
---|
253 | const SaturationsPositionsPerNeighbor_t &saturationpositions =
|
---|
254 | mapiter->second;
|
---|
255 |
|
---|
256 | // go through all cut bonds for this atom
|
---|
257 | for (BondList::const_iterator bonditer = atomiter->second.begin();
|
---|
258 | bonditer != atomiter->second.end(); ++bonditer) {
|
---|
259 | atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
|
---|
260 |
|
---|
261 | // get positions from global map
|
---|
262 | SaturationsPositionsPerNeighbor_t::const_iterator positionsiter =
|
---|
263 | saturationpositions.find(OtherWalker->getId());
|
---|
264 | ASSERT(positionsiter != saturationpositions.end(),
|
---|
265 | "SaturatedFragment::saturate() - no information on bond neighbor "
|
---|
266 | +toString(*OtherWalker)+" to atom "+toString(*Walker));
|
---|
267 | ASSERT(!positionsiter->second.empty(),
|
---|
268 | "SaturatedFragment::saturate() - no positions for saturating bond to"
|
---|
269 | +toString(*OtherWalker)+" to atom "+toString(*Walker));
|
---|
270 |
|
---|
271 | // get typical bond distance from elements database
|
---|
272 | double BondDistance = Walker->getType()->getHBondDistance(positionsiter->second.size()-1);
|
---|
273 | if (BondDistance < 0.) {
|
---|
274 | ELOG(2, "saturateAtoms() - no typical hydrogen bond distance of degree "
|
---|
275 | +toString(positionsiter->second.size())+" for element "
|
---|
276 | +toString(Walker->getType()->getName()));
|
---|
277 | // try bond degree 1 distance
|
---|
278 | BondDistance = Walker->getType()->getHBondDistance(1-1);
|
---|
279 | if (BondDistance < 0.) {
|
---|
280 | ELOG(1, "saturateAtoms() - no typical hydrogen bond distance for element "
|
---|
281 | +toString(Walker->getType()->getName()));
|
---|
282 | BondDistance = 1.;
|
---|
283 | }
|
---|
284 | }
|
---|
285 | ASSERT( BondDistance > 0.,
|
---|
286 | "SaturatedFragment::saturate() - negative bond distance");
|
---|
287 |
|
---|
288 | // place hydrogen at each point
|
---|
289 | LOG(4, "DEBUG: Places to saturate for atom " << *OtherWalker
|
---|
290 | << " are " << positionsiter->second);
|
---|
291 | atom * const father = Walker;
|
---|
292 | for (SaturationsPositions_t::const_iterator positer = positionsiter->second.begin();
|
---|
293 | positer != positionsiter->second.end(); ++positer) {
|
---|
294 | const atom& hydrogen =
|
---|
295 | setHydrogenReplacement(Walker, *positer, BondDistance, father);
|
---|
296 | FullMolecule.insert(hydrogen.getId());
|
---|
297 | }
|
---|
298 | }
|
---|
299 | }
|
---|
300 | } else
|
---|
301 | LOG(3, "INFO: We are not saturating cut bonds.");
|
---|
302 | }
|
---|
303 |
|
---|
304 | const atom& SaturatedFragment::setHydrogenReplacement(
|
---|
305 | const atom * const _OwnerAtom,
|
---|
306 | const Vector &_position,
|
---|
307 | const double _distance,
|
---|
308 | atom * const _father)
|
---|
309 | {
|
---|
310 | atom * const _atom = hydrogens.leaseHydrogen(); // new atom
|
---|
311 | _atom->setPosition( _OwnerAtom->getPosition() + _distance * _position );
|
---|
312 | // always set as fixed ion (not moving during molecular dynamics simulation)
|
---|
313 | _atom->setFixedIon(true);
|
---|
314 | // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
|
---|
315 | _atom->father = _father;
|
---|
316 | SaturationHydrogens.insert(_atom->getId());
|
---|
317 |
|
---|
318 | return *_atom;
|
---|
319 | }
|
---|
320 |
|
---|
321 | bool SaturatedFragment::saturateAtom(
|
---|
322 | atom * const _atom,
|
---|
323 | const BondList &_cutbonds)
|
---|
324 | {
|
---|
325 | // go through each bond and replace
|
---|
326 | for (BondList::const_iterator bonditer = _cutbonds.begin();
|
---|
327 | bonditer != _cutbonds.end(); ++bonditer) {
|
---|
328 | atom * const OtherWalker = (*bonditer)->GetOtherAtom(_atom);
|
---|
329 | if (!AddHydrogenReplacementAtom(
|
---|
330 | (*bonditer),
|
---|
331 | _atom,
|
---|
332 | OtherWalker,
|
---|
333 | World::getInstance().getConfig()->IsAngstroem == 1))
|
---|
334 | return false;
|
---|
335 | }
|
---|
336 | return true;
|
---|
337 | }
|
---|
338 |
|
---|
339 | bool SaturatedFragment::OutputConfig(
|
---|
340 | std::ostream &out,
|
---|
341 | const ParserTypes _type) const
|
---|
342 | {
|
---|
343 | // gather all atoms in a vector
|
---|
344 | std::vector<atom *> atoms;
|
---|
345 | for (KeySet::const_iterator iter = FullMolecule.begin();
|
---|
346 | iter != FullMolecule.end();
|
---|
347 | ++iter) {
|
---|
348 | atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
|
---|
349 | ASSERT( Walker != NULL,
|
---|
350 | "SaturatedFragment::OutputConfig() - id "
|
---|
351 | +toString(*iter)+" is unknown to World.");
|
---|
352 | atoms.push_back(Walker);
|
---|
353 | }
|
---|
354 |
|
---|
355 | // TODO: ScanForPeriodicCorrection() is missing so far!
|
---|
356 | // note however that this is not straight-forward for the following reasons:
|
---|
357 | // - we do not copy all atoms anymore, hence we are forced to shift the real
|
---|
358 | // atoms hither and back again
|
---|
359 | // - we use a long-range potential that supports periodic boundary conditions.
|
---|
360 | // Hence, there we would like the original configuration (split through the
|
---|
361 | // the periodic boundaries). Otherwise, we would have to shift (and probably
|
---|
362 | // interpolate) the potential with OBCs applying.
|
---|
363 |
|
---|
364 | // list atoms in fragment for debugging
|
---|
365 | {
|
---|
366 | std::stringstream output;
|
---|
367 | output << "INFO: Contained atoms: ";
|
---|
368 | for (std::vector<atom *>::const_iterator iter = atoms.begin();
|
---|
369 | iter != atoms.end(); ++iter) {
|
---|
370 | output << (*iter)->getName() << " ";
|
---|
371 | }
|
---|
372 | LOG(3, output.str());
|
---|
373 | }
|
---|
374 |
|
---|
375 | // store to stream via FragmentParser
|
---|
376 | const bool intermediateResult =
|
---|
377 | FormatParserStorage::getInstance().save(
|
---|
378 | out,
|
---|
379 | FormatParserStorage::getInstance().getSuffixFromType(_type),
|
---|
380 | atoms);
|
---|
381 |
|
---|
382 | return intermediateResult;
|
---|
383 | }
|
---|
384 |
|
---|
385 | atom * const SaturatedFragment::getHydrogenReplacement(atom * const replacement)
|
---|
386 | {
|
---|
387 | atom * const _atom = hydrogens.leaseHydrogen(); // new atom
|
---|
388 | _atom->setAtomicVelocity(replacement->getAtomicVelocity()); // copy velocity
|
---|
389 | _atom->setFixedIon(replacement->getFixedIon());
|
---|
390 | // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
|
---|
391 | _atom->father = replacement;
|
---|
392 | SaturationHydrogens.insert(_atom->getId());
|
---|
393 | return _atom;
|
---|
394 | }
|
---|
395 |
|
---|
396 | bool SaturatedFragment::AddHydrogenReplacementAtom(
|
---|
397 | bond::ptr TopBond,
|
---|
398 | atom *Origin,
|
---|
399 | atom *Replacement,
|
---|
400 | bool IsAngstroem)
|
---|
401 | {
|
---|
402 | // Info info(__func__);
|
---|
403 | bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
|
---|
404 | double bondlength; // bond length of the bond to be replaced/cut
|
---|
405 | double bondangle; // bond angle of the bond to be replaced/cut
|
---|
406 | double BondRescale; // rescale value for the hydrogen bond length
|
---|
407 | bond::ptr FirstBond;
|
---|
408 | bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
|
---|
409 | atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
|
---|
410 | double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
|
---|
411 | Vector Orthovector1, Orthovector2; // temporary vectors in coordination construction
|
---|
412 | Vector InBondvector; // vector in direction of *Bond
|
---|
413 | const RealSpaceMatrix &matrix = World::getInstance().getDomain().getM();
|
---|
414 | bond::ptr Binder;
|
---|
415 |
|
---|
416 | // create vector in direction of bond
|
---|
417 | InBondvector = Replacement->getPosition() - Origin->getPosition();
|
---|
418 | bondlength = InBondvector.Norm();
|
---|
419 |
|
---|
420 | // is greater than typical bond distance? Then we have to correct periodically
|
---|
421 | // the problem is not the H being out of the box, but InBondvector have the wrong direction
|
---|
422 | // due to Replacement or Origin being on the wrong side!
|
---|
423 | const BondGraph * const BG = World::getInstance().getBondGraph();
|
---|
424 | const range<double> MinMaxBondDistance(
|
---|
425 | BG->getMinMaxDistance(Origin,Replacement));
|
---|
426 | if (!MinMaxBondDistance.isInRange(bondlength)) {
|
---|
427 | // LOG(4, "InBondvector is: " << InBondvector << ".");
|
---|
428 | Orthovector1.Zero();
|
---|
429 | for (int i=NDIM;i--;) {
|
---|
430 | l = Replacement->at(i) - Origin->at(i);
|
---|
431 | if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
|
---|
432 | Orthovector1[i] = (l < 0) ? -1. : +1.;
|
---|
433 | } // (signs are correct, was tested!)
|
---|
434 | }
|
---|
435 | Orthovector1 *= matrix;
|
---|
436 | InBondvector -= Orthovector1; // subtract just the additional translation
|
---|
437 | bondlength = InBondvector.Norm();
|
---|
438 | // LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
|
---|
439 | } // periodic correction finished
|
---|
440 |
|
---|
441 | InBondvector.Normalize();
|
---|
442 | // get typical bond length and store as scale factor for later
|
---|
443 | ASSERT(Origin->getType() != NULL,
|
---|
444 | "SaturatedFragment::AddHydrogenReplacementAtom() - element of Origin is not given.");
|
---|
445 | BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree()-1);
|
---|
446 | if (BondRescale == -1) {
|
---|
447 | ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
|
---|
448 | BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
|
---|
449 | if (BondRescale == -1) {
|
---|
450 | ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
|
---|
451 | return false;
|
---|
452 | BondRescale = bondlength;
|
---|
453 | }
|
---|
454 | } else {
|
---|
455 | if (!IsAngstroem)
|
---|
456 | BondRescale /= (1.*AtomicLengthToAngstroem);
|
---|
457 | }
|
---|
458 |
|
---|
459 | // discern single, double and triple bonds
|
---|
460 | switch(TopBond->getDegree()) {
|
---|
461 | case 1:
|
---|
462 | // check whether replacement has been an excluded hydrogen
|
---|
463 | if (Replacement->getType()->getAtomicNumber() == HydrogenPool::HYDROGEN) { // neither rescale nor replace if it's already hydrogen
|
---|
464 | FirstOtherAtom = Replacement;
|
---|
465 | BondRescale = bondlength;
|
---|
466 | } else {
|
---|
467 | FirstOtherAtom = getHydrogenReplacement(Replacement);
|
---|
468 | InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length
|
---|
469 | FirstOtherAtom->setPosition(Origin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
|
---|
470 | }
|
---|
471 | FullMolecule.insert(FirstOtherAtom->getId());
|
---|
472 | // LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
|
---|
473 | break;
|
---|
474 | case 2:
|
---|
475 | {
|
---|
476 | // determine two other bonds (warning if there are more than two other) plus valence sanity check
|
---|
477 | const BondList& ListOfBonds = Origin->getListOfBonds();
|
---|
478 | for (BondList::const_iterator Runner = ListOfBonds.begin();
|
---|
479 | Runner != ListOfBonds.end();
|
---|
480 | ++Runner) {
|
---|
481 | if ((*Runner) != TopBond) {
|
---|
482 | if (FirstBond == NULL) {
|
---|
483 | FirstBond = (*Runner);
|
---|
484 | FirstOtherAtom = (*Runner)->GetOtherAtom(Origin);
|
---|
485 | } else if (SecondBond == NULL) {
|
---|
486 | SecondBond = (*Runner);
|
---|
487 | SecondOtherAtom = (*Runner)->GetOtherAtom(Origin);
|
---|
488 | } else {
|
---|
489 | ELOG(2, "Detected more than four bonds for atom " << Origin->getName());
|
---|
490 | }
|
---|
491 | }
|
---|
492 | }
|
---|
493 | }
|
---|
494 | 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)
|
---|
495 | SecondBond = TopBond;
|
---|
496 | SecondOtherAtom = Replacement;
|
---|
497 | }
|
---|
498 | if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
|
---|
499 | // 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.");
|
---|
500 |
|
---|
501 | // determine the plane of these two with the *origin
|
---|
502 | try {
|
---|
503 | Orthovector1 = Plane(Origin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
|
---|
504 | }
|
---|
505 | catch(LinearDependenceException &excp){
|
---|
506 | LOG(0, boost::diagnostic_information(excp));
|
---|
507 | // TODO: figure out what to do with the Orthovector in this case
|
---|
508 | AllWentWell = false;
|
---|
509 | }
|
---|
510 | } else {
|
---|
511 | Orthovector1.GetOneNormalVector(InBondvector);
|
---|
512 | }
|
---|
513 | //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
|
---|
514 | // orthogonal vector and bond vector between origin and replacement form the new plane
|
---|
515 | Orthovector1.MakeNormalTo(InBondvector);
|
---|
516 | Orthovector1.Normalize();
|
---|
517 | //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
|
---|
518 |
|
---|
519 | // create the two Hydrogens ...
|
---|
520 | FirstOtherAtom = getHydrogenReplacement(Replacement);
|
---|
521 | SecondOtherAtom = getHydrogenReplacement(Replacement);
|
---|
522 | FullMolecule.insert(FirstOtherAtom->getId());
|
---|
523 | FullMolecule.insert(SecondOtherAtom->getId());
|
---|
524 | bondangle = Origin->getType()->getHBondAngle(1);
|
---|
525 | if (bondangle == -1) {
|
---|
526 | ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
|
---|
527 | return false;
|
---|
528 | bondangle = 0;
|
---|
529 | }
|
---|
530 | bondangle *= M_PI/180./2.;
|
---|
531 | // LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
|
---|
532 | // LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
|
---|
533 | FirstOtherAtom->Zero();
|
---|
534 | SecondOtherAtom->Zero();
|
---|
535 | for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
|
---|
536 | FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
|
---|
537 | SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
|
---|
538 | }
|
---|
539 | FirstOtherAtom->Scale(BondRescale); // rescale by correct BondDistance
|
---|
540 | SecondOtherAtom->Scale(BondRescale);
|
---|
541 | //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
|
---|
542 | *FirstOtherAtom += Origin->getPosition();
|
---|
543 | *SecondOtherAtom += Origin->getPosition();
|
---|
544 | // ... and add to molecule
|
---|
545 | // LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
|
---|
546 | // LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
|
---|
547 | break;
|
---|
548 | case 3:
|
---|
549 | // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
|
---|
550 | FirstOtherAtom = getHydrogenReplacement(Replacement);
|
---|
551 | SecondOtherAtom = getHydrogenReplacement(Replacement);
|
---|
552 | ThirdOtherAtom = getHydrogenReplacement(Replacement);
|
---|
553 | FullMolecule.insert(FirstOtherAtom->getId());
|
---|
554 | FullMolecule.insert(SecondOtherAtom->getId());
|
---|
555 | FullMolecule.insert(ThirdOtherAtom->getId());
|
---|
556 |
|
---|
557 | // we need to vectors orthonormal the InBondvector
|
---|
558 | AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
|
---|
559 | // LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
|
---|
560 | try{
|
---|
561 | Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
|
---|
562 | }
|
---|
563 | catch(LinearDependenceException &excp) {
|
---|
564 | LOG(0, boost::diagnostic_information(excp));
|
---|
565 | AllWentWell = false;
|
---|
566 | }
|
---|
567 | // LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
|
---|
568 |
|
---|
569 | // create correct coordination for the three atoms
|
---|
570 | alpha = (Origin->getType()->getHBondAngle(2))/180.*M_PI/2.; // retrieve triple bond angle from database
|
---|
571 | l = BondRescale; // desired bond length
|
---|
572 | b = 2.*l*sin(alpha); // base length of isosceles triangle
|
---|
573 | d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
|
---|
574 | f = b/sqrt(3.); // length for Orthvector1
|
---|
575 | g = b/2.; // length for Orthvector2
|
---|
576 | // LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
|
---|
577 | // 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));
|
---|
578 | factors[0] = d;
|
---|
579 | factors[1] = f;
|
---|
580 | factors[2] = 0.;
|
---|
581 | FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
|
---|
582 | factors[1] = -0.5*f;
|
---|
583 | factors[2] = g;
|
---|
584 | SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
|
---|
585 | factors[2] = -g;
|
---|
586 | ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
|
---|
587 |
|
---|
588 | // rescale each to correct BondDistance
|
---|
589 | // FirstOtherAtom->x.Scale(&BondRescale);
|
---|
590 | // SecondOtherAtom->x.Scale(&BondRescale);
|
---|
591 | // ThirdOtherAtom->x.Scale(&BondRescale);
|
---|
592 |
|
---|
593 | // and relative to *origin atom
|
---|
594 | *FirstOtherAtom += Origin->getPosition();
|
---|
595 | *SecondOtherAtom += Origin->getPosition();
|
---|
596 | *ThirdOtherAtom += Origin->getPosition();
|
---|
597 |
|
---|
598 | // ... and add to molecule
|
---|
599 | // LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
|
---|
600 | // LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
|
---|
601 | // LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
|
---|
602 | break;
|
---|
603 | default:
|
---|
604 | ELOG(1, "BondDegree does not state single, double or triple bond!");
|
---|
605 | AllWentWell = false;
|
---|
606 | break;
|
---|
607 | }
|
---|
608 |
|
---|
609 | return AllWentWell;
|
---|
610 | };
|
---|