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

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 9e3fca was 9d3264, checked in by Frederik Heber <heber@…>, 11 years ago

Refactored SaturatedFragment::saturate() using list of cut bonds.

  • Property mode set to 100644
File size: 19.1 KB
RevLine 
[7d5fcd]1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2013 University of Bonn. All rights reserved.
[5aaa43]5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
[7d5fcd]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
[c39675]40#include <cmath>
41#include <iostream>
42
[7d5fcd]43#include "CodePatterns/Assert.hpp"
[c39675]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"
[7d5fcd]50
[c39675]51#include "Atom/atom.hpp"
52#include "Bond/bond.hpp"
53#include "config.hpp"
54#include "Descriptors/AtomIdDescriptor.hpp"
[7d5fcd]55#include "Fragmentation/Exporters/HydrogenPool.hpp"
[c39675]56#include "Fragmentation/HydrogenSaturation_enum.hpp"
57#include "Graph/BondGraph.hpp"
58#include "World.hpp"
[7d5fcd]59
60SaturatedFragment::SaturatedFragment(
61 const KeySet &_set,
62 KeySetsInUse_t &_container,
[c39675]63 HydrogenPool &_hydrogens,
64 const enum HydrogenTreatment _treatment,
65 const enum HydrogenSaturation _saturation) :
[7d5fcd]66 container(_container),
67 set(_set),
68 hydrogens(_hydrogens),
[c39675]69 FullMolecule(set),
70 treatment(_treatment),
71 saturation(_saturation)
[7d5fcd]72{
73 // add to in-use container
74 ASSERT( container.find(set) == container.end(),
75 "SaturatedFragment::SaturatedFragment() - the set "
76 +toString(set)+" is already marked as in use.");
77 container.insert(set);
78
[c39675]79 // prepare saturation hydrogens
80 saturate();
[7d5fcd]81}
82
83SaturatedFragment::~SaturatedFragment()
84{
85 // release all saturation hydrogens if present
86 for (KeySet::iterator iter = SaturationHydrogens.begin();
87 !SaturationHydrogens.empty();
88 iter = SaturationHydrogens.begin()) {
89 hydrogens.releaseHydrogen(*iter);
90 SaturationHydrogens.erase(iter);
91 }
92
93 // remove ourselves from in-use container
94 KeySetsInUse_t::iterator iter = container.find(set);
95 ASSERT( container.find(set) != container.end(),
96 "SaturatedFragment::SaturatedFragment() - the set "
97 +toString(set)+" is not marked as in use.");
98 container.erase(iter);
99}
[c39675]100
101void SaturatedFragment::saturate()
102{
[9d3264]103 // so far, we just have a set of keys. Hence, convert to atom refs
104 // and gather all atoms in a vector
[c39675]105 std::vector<atom *> atoms;
106 for (KeySet::const_iterator iter = FullMolecule.begin();
107 iter != FullMolecule.end();
108 ++iter) {
109 atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
110 ASSERT( Walker != NULL,
111 "SaturatedFragment::OutputConfig() - id "
112 +toString(*iter)+" is unknown to World.");
113 atoms.push_back(Walker);
114 }
115
116// bool LonelyFlag = false;
[9d3264]117 // go through each atom of the fragment and gather all cut bonds in list
118 typedef std::map<atom *, BondList > CutBonds_t;
119 CutBonds_t CutBonds;
[c39675]120 for (std::vector<atom *>::const_iterator iter = atoms.begin();
121 iter != atoms.end();
122 ++iter) {
123 atom * const Walker = *iter;
124
125 // go through all bonds
126 const BondList& ListOfBonds = Walker->getListOfBonds();
127 for (BondList::const_iterator BondRunner = ListOfBonds.begin();
128 BondRunner != ListOfBonds.end();
129 ++BondRunner) {
130 atom * const OtherWalker = (*BondRunner)->GetOtherAtom(Walker);
[9d3264]131 // if other atom is in key set
[c39675]132 if (set.find(OtherWalker->getId()) != set.end()) {
133 LOG(4, "DEBUG: Walker " << *Walker << " is bound to " << *OtherWalker << ".");
134// if (OtherWalker->getId() > Walker->getId()) { // add bond (Nr check is for adding only one of both variants: ab, ba)
135//// std::stringstream output;
136//// output << "ACCEPT: Adding Bond: "
[1f693d]137// output << Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
[c39675]138//// LOG(3, output.str());
139// //NumBonds[(*iter)->getNr()]++;
140// } else {
141//// LOG(3, "REJECY: Not adding bond, labels in wrong order.");
142// }
143// LonelyFlag = false;
144 } else {
145 LOG(4, "DEBUG: Walker " << *Walker << " is bound to "
146 << *OtherWalker << ", who is not in this fragment molecule.");
147 if (saturation == DoSaturate) {
148// LOG(3, "ACCEPT: Adding Hydrogen to " << (*iter)->Name << " and a bond in between.");
[9d3264]149 if (CutBonds.count(Walker) == 0)
150 CutBonds.insert( std::make_pair(Walker, BondList() ));
151 CutBonds[Walker].push_back(*BondRunner);
[c39675]152 }
153// } else if ((treatment == ExcludeHydrogen) && (OtherWalker->getElementNo() == (atomicNumber_t)1)) {
154// // just copy the atom if it's a hydrogen
155// atom * const OtherWalker = Leaf->AddCopyAtom(OtherWalker);
[1f693d]156// Leaf->AddBond((*iter), OtherWalker, (*BondRunner)->getDegree());
[c39675]157// }
[1f693d]158 //NumBonds[(*iter)->getNr()] += Binder->getDegree();
[c39675]159 }
160 }
161 }
[9d3264]162
163 // go through all cut bonds and replace with a hydrogen
164 for (CutBonds_t::const_iterator atomiter = CutBonds.begin();
165 atomiter != CutBonds.end(); ++atomiter) {
166 atom * const Walker = atomiter->first;
167 // go through each bond and replace
168 for (BondList::const_iterator bonditer = atomiter->second.begin();
169 bonditer != atomiter->second.end(); ++bonditer) {
170 atom * const OtherWalker = (*bonditer)->GetOtherAtom(Walker);
171 if (!AddHydrogenReplacementAtom(
172 (*bonditer),
173 Walker,
174 OtherWalker,
175 World::getInstance().getConfig()->IsAngstroem == 1))
176 exit(1);
177 }
178 }
[c39675]179}
180
181bool SaturatedFragment::OutputConfig(
182 std::ostream &out,
183 const ParserTypes _type) const
184{
185 // gather all atoms in a vector
186 std::vector<atom *> atoms;
187 for (KeySet::const_iterator iter = FullMolecule.begin();
188 iter != FullMolecule.end();
189 ++iter) {
190 atom * const Walker = World::getInstance().getAtom(AtomById(*iter));
191 ASSERT( Walker != NULL,
192 "SaturatedFragment::OutputConfig() - id "
193 +toString(*iter)+" is unknown to World.");
194 atoms.push_back(Walker);
195 }
196
197 // TODO: ScanForPeriodicCorrection() is missing so far!
198 // note however that this is not straight-forward for the following reasons:
199 // - we do not copy all atoms anymore, hence we are forced to shift the real
200 // atoms hither and back again
201 // - we use a long-range potential that supports periodic boundary conditions.
202 // Hence, there we would like the original configuration (split through the
203 // the periodic boundaries). Otherwise, we would have to shift (and probably
204 // interpolate) the potential with OBCs applying.
205
206 // list atoms in fragment for debugging
207 {
208 std::stringstream output;
209 output << "INFO: Contained atoms: ";
210 for (std::vector<atom *>::const_iterator iter = atoms.begin();
211 iter != atoms.end(); ++iter) {
212 output << (*iter)->getName() << " ";
213 }
214 LOG(3, output.str());
215 }
216
217 // store to stream via FragmentParser
218 const bool intermediateResult =
219 FormatParserStorage::getInstance().save(
220 out,
221 FormatParserStorage::getInstance().getSuffixFromType(_type),
222 atoms);
223
224 return intermediateResult;
225}
226
227atom * const SaturatedFragment::getHydrogenReplacement(atom * const replacement)
228{
229 atom * const _atom = hydrogens.leaseHydrogen(); // new atom
230 _atom->setAtomicVelocity(replacement->getAtomicVelocity()); // copy velocity
231 _atom->setFixedIon(replacement->getFixedIon());
232 // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
233 _atom->father = replacement;
234 SaturationHydrogens.insert(_atom->getId());
235 return _atom;
236}
237
238bool SaturatedFragment::AddHydrogenReplacementAtom(
239 bond::ptr TopBond,
240 atom *Origin,
241 atom *Replacement,
242 bool IsAngstroem)
243{
244// Info info(__func__);
245 bool AllWentWell = true; // flag gathering the boolean return value of molecule::AddAtom and other functions, as return value on exit
246 double bondlength; // bond length of the bond to be replaced/cut
247 double bondangle; // bond angle of the bond to be replaced/cut
248 double BondRescale; // rescale value for the hydrogen bond length
249 bond::ptr FirstBond;
250 bond::ptr SecondBond; // Other bonds in double bond case to determine "other" plane
251 atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
252 double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination
253 Vector Orthovector1, Orthovector2; // temporary vectors in coordination construction
254 Vector InBondvector; // vector in direction of *Bond
255 const RealSpaceMatrix &matrix = World::getInstance().getDomain().getM();
256 bond::ptr Binder;
257
258 // create vector in direction of bond
259 InBondvector = Replacement->getPosition() - Origin->getPosition();
260 bondlength = InBondvector.Norm();
261
262 // is greater than typical bond distance? Then we have to correct periodically
263 // the problem is not the H being out of the box, but InBondvector have the wrong direction
264 // due to Replacement or Origin being on the wrong side!
265 const BondGraph * const BG = World::getInstance().getBondGraph();
266 const range<double> MinMaxBondDistance(
267 BG->getMinMaxDistance(Origin,Replacement));
268 if (!MinMaxBondDistance.isInRange(bondlength)) {
269// LOG(4, "InBondvector is: " << InBondvector << ".");
270 Orthovector1.Zero();
271 for (int i=NDIM;i--;) {
272 l = Replacement->at(i) - Origin->at(i);
273 if (fabs(l) > MinMaxBondDistance.last) { // is component greater than bond distance (check against min not useful here)
274 Orthovector1[i] = (l < 0) ? -1. : +1.;
275 } // (signs are correct, was tested!)
276 }
277 Orthovector1 *= matrix;
278 InBondvector -= Orthovector1; // subtract just the additional translation
279 bondlength = InBondvector.Norm();
280// LOG(4, "INFO: Corrected InBondvector is now: " << InBondvector << ".");
281 } // periodic correction finished
282
283 InBondvector.Normalize();
284 // get typical bond length and store as scale factor for later
285 ASSERT(Origin->getType() != NULL,
286 "SaturatedFragment::AddHydrogenReplacementAtom() - element of Origin is not given.");
[1f693d]287 BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree()-1);
[c39675]288 if (BondRescale == -1) {
[1f693d]289 ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
[3fbdca]290 BondRescale = Origin->getType()->getHBondDistance(TopBond->getDegree());
291 if (BondRescale == -1) {
292 ELOG(1, "There is no typical hydrogen bond distance in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of any degree!");
293 return false;
294 BondRescale = bondlength;
295 }
[c39675]296 } else {
297 if (!IsAngstroem)
298 BondRescale /= (1.*AtomicLengthToAngstroem);
299 }
300
301 // discern single, double and triple bonds
[1f693d]302 switch(TopBond->getDegree()) {
[c39675]303 case 1:
304 // check whether replacement has been an excluded hydrogen
305 if (Replacement->getType()->getAtomicNumber() == HydrogenPool::HYDROGEN) { // neither rescale nor replace if it's already hydrogen
306 FirstOtherAtom = Replacement;
307 BondRescale = bondlength;
308 } else {
309 FirstOtherAtom = getHydrogenReplacement(Replacement);
310 InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length
311 FirstOtherAtom->setPosition(Origin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
312 }
313 FullMolecule.insert(FirstOtherAtom->getId());
314// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
315 break;
316 case 2:
317 {
318 // determine two other bonds (warning if there are more than two other) plus valence sanity check
319 const BondList& ListOfBonds = Origin->getListOfBonds();
320 for (BondList::const_iterator Runner = ListOfBonds.begin();
321 Runner != ListOfBonds.end();
322 ++Runner) {
323 if ((*Runner) != TopBond) {
324 if (FirstBond == NULL) {
325 FirstBond = (*Runner);
326 FirstOtherAtom = (*Runner)->GetOtherAtom(Origin);
327 } else if (SecondBond == NULL) {
328 SecondBond = (*Runner);
329 SecondOtherAtom = (*Runner)->GetOtherAtom(Origin);
330 } else {
331 ELOG(2, "Detected more than four bonds for atom " << Origin->getName());
332 }
333 }
334 }
335 }
336 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)
337 SecondBond = TopBond;
338 SecondOtherAtom = Replacement;
339 }
340 if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
341// 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.");
342
343 // determine the plane of these two with the *origin
344 try {
345 Orthovector1 = Plane(Origin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
346 }
347 catch(LinearDependenceException &excp){
348 LOG(0, boost::diagnostic_information(excp));
349 // TODO: figure out what to do with the Orthovector in this case
350 AllWentWell = false;
351 }
352 } else {
353 Orthovector1.GetOneNormalVector(InBondvector);
354 }
355 //LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
356 // orthogonal vector and bond vector between origin and replacement form the new plane
357 Orthovector1.MakeNormalTo(InBondvector);
358 Orthovector1.Normalize();
359 //LOG(3, "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << ".");
360
361 // create the two Hydrogens ...
362 FirstOtherAtom = getHydrogenReplacement(Replacement);
363 SecondOtherAtom = getHydrogenReplacement(Replacement);
364 FullMolecule.insert(FirstOtherAtom->getId());
365 FullMolecule.insert(SecondOtherAtom->getId());
366 bondangle = Origin->getType()->getHBondAngle(1);
367 if (bondangle == -1) {
[1f693d]368 ELOG(1, "There is no typical hydrogen bond angle in replacing bond (" << Origin->getName() << "<->" << Replacement->getName() << ") of degree " << TopBond->getDegree() << "!");
[c39675]369 return false;
370 bondangle = 0;
371 }
372 bondangle *= M_PI/180./2.;
373// LOG(3, "INFO: ReScaleCheck: InBondvector " << InBondvector << ", " << Orthovector1 << ".");
374// LOG(3, "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle));
375 FirstOtherAtom->Zero();
376 SecondOtherAtom->Zero();
377 for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
378 FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
379 SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
380 }
381 FirstOtherAtom->Scale(BondRescale); // rescale by correct BondDistance
382 SecondOtherAtom->Scale(BondRescale);
383 //LOG(3, "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << ".");
384 *FirstOtherAtom += Origin->getPosition();
385 *SecondOtherAtom += Origin->getPosition();
386 // ... and add to molecule
387// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
388// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
389 break;
390 case 3:
391 // take the "usual" tetraoidal angle and add the three Hydrogen in direction of the bond (height of the tetraoid)
392 FirstOtherAtom = getHydrogenReplacement(Replacement);
393 SecondOtherAtom = getHydrogenReplacement(Replacement);
394 ThirdOtherAtom = getHydrogenReplacement(Replacement);
395 FullMolecule.insert(FirstOtherAtom->getId());
396 FullMolecule.insert(SecondOtherAtom->getId());
397 FullMolecule.insert(ThirdOtherAtom->getId());
398
399 // we need to vectors orthonormal the InBondvector
400 AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
401// LOG(3, "INFO: Orthovector1: " << Orthovector1 << ".");
402 try{
403 Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
404 }
405 catch(LinearDependenceException &excp) {
406 LOG(0, boost::diagnostic_information(excp));
407 AllWentWell = false;
408 }
409// LOG(3, "INFO: Orthovector2: " << Orthovector2 << ".")
410
411 // create correct coordination for the three atoms
412 alpha = (Origin->getType()->getHBondAngle(2))/180.*M_PI/2.; // retrieve triple bond angle from database
413 l = BondRescale; // desired bond length
414 b = 2.*l*sin(alpha); // base length of isosceles triangle
415 d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector
416 f = b/sqrt(3.); // length for Orthvector1
417 g = b/2.; // length for Orthvector2
418// LOG(3, "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", ");
419// 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));
420 factors[0] = d;
421 factors[1] = f;
422 factors[2] = 0.;
423 FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
424 factors[1] = -0.5*f;
425 factors[2] = g;
426 SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
427 factors[2] = -g;
428 ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
429
430 // rescale each to correct BondDistance
431// FirstOtherAtom->x.Scale(&BondRescale);
432// SecondOtherAtom->x.Scale(&BondRescale);
433// ThirdOtherAtom->x.Scale(&BondRescale);
434
435 // and relative to *origin atom
436 *FirstOtherAtom += Origin->getPosition();
437 *SecondOtherAtom += Origin->getPosition();
438 *ThirdOtherAtom += Origin->getPosition();
439
440 // ... and add to molecule
441// LOG(4, "INFO: Added " << *FirstOtherAtom << " at: " << FirstOtherAtom->x << ".");
442// LOG(4, "INFO: Added " << *SecondOtherAtom << " at: " << SecondOtherAtom->x << ".");
443// LOG(4, "INFO: Added " << *ThirdOtherAtom << " at: " << ThirdOtherAtom->x << ".");
444 break;
445 default:
446 ELOG(1, "BondDegree does not state single, double or triple bond!");
447 AllWentWell = false;
448 break;
449 }
450
451 return AllWentWell;
452};
Note: See TracBrowser for help on using the repository browser.