Changes in src/molecule.cpp [d74077:952f38]
- File:
-
- 1 edited
-
src/molecule.cpp (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
rd74077 r952f38 24 24 #include "element.hpp" 25 25 #include "graph.hpp" 26 #include " helpers.hpp"26 #include "Helpers/helpers.hpp" 27 27 #include "leastsquaremin.hpp" 28 28 #include "linkedcell.hpp" 29 29 #include "lists.hpp" 30 #include " log.hpp"30 #include "Helpers/Log.hpp" 31 31 #include "molecule.hpp" 32 32 … … 34 34 #include "stackclass.hpp" 35 35 #include "tesselation.hpp" 36 #include " vector.hpp"37 #include " Matrix.hpp"36 #include "LinearAlgebra/Vector.hpp" 37 #include "LinearAlgebra/Matrix.hpp" 38 38 #include "World.hpp" 39 39 #include "Box.hpp" 40 #include " Plane.hpp"40 #include "LinearAlgebra/Plane.hpp" 41 41 #include "Exceptions/LinearDependenceException.hpp" 42 42 … … 206 206 if (pointer != NULL) { 207 207 pointer->sort = &pointer->nr; 208 if (pointer-> getType()!= NULL) {209 formula += pointer-> getType();210 if (pointer-> getType()->Z != 1)208 if (pointer->type != NULL) { 209 formula += pointer->type; 210 if (pointer->type->Z != 1) 211 211 NoNonHydrogen++; 212 212 if(pointer->getName() == "Unknown"){ 213 213 stringstream sstr; 214 sstr << pointer-> getType()->symbol << pointer->nr+1;214 sstr << pointer->type->symbol << pointer->nr+1; 215 215 pointer->setName(sstr.str()); 216 216 } … … 236 236 walker->nr = last_atom++; // increase number within molecule 237 237 insert(walker); 238 if ((pointer-> getType() != NULL) && (pointer->getType()->Z != 1))238 if ((pointer->type != NULL) && (pointer->type->Z != 1)) 239 239 NoNonHydrogen++; 240 240 retval=walker; … … 294 294 // Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl; 295 295 // create vector in direction of bond 296 InBondvector = TopReplacement-> getPosition() - TopOrigin->getPosition();296 InBondvector = TopReplacement->x - TopOrigin->x; 297 297 bondlength = InBondvector.Norm(); 298 298 … … 306 306 Orthovector1.Zero(); 307 307 for (int i=NDIM;i--;) { 308 l = TopReplacement-> at(i) - TopOrigin->at(i);308 l = TopReplacement->x[i] - TopOrigin->x[i]; 309 309 if (fabs(l) > BondDistance) { // is component greater than bond distance 310 310 Orthovector1[i] = (l < 0) ? -1. : +1.; … … 321 321 InBondvector.Normalize(); 322 322 // get typical bond length and store as scale factor for later 323 ASSERT(TopOrigin-> getType()!= NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");324 BondRescale = TopOrigin-> getType()->HBondDistance[TopBond->BondDegree-1];323 ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given."); 324 BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1]; 325 325 if (BondRescale == -1) { 326 326 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl); … … 336 336 case 1: 337 337 FirstOtherAtom = World::getInstance().createAtom(); // new atom 338 FirstOtherAtom-> setType(1); // element is Hydrogen339 FirstOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity338 FirstOtherAtom->type = elemente->FindElement(1); // element is Hydrogen 339 FirstOtherAtom->v = TopReplacement->v; // copy velocity 340 340 FirstOtherAtom->FixedIon = TopReplacement->FixedIon; 341 if (TopReplacement-> getType()->Z == 1) { // neither rescale nor replace if it's already hydrogen341 if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen 342 342 FirstOtherAtom->father = TopReplacement; 343 343 BondRescale = bondlength; … … 346 346 } 347 347 InBondvector *= BondRescale; // rescale the distance vector to Hydrogen bond length 348 FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom 348 FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ... 349 FirstOtherAtom->x += InBondvector; // ... and add distance vector to replacement atom 349 350 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom); 350 351 // Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: "; … … 379 380 // determine the plane of these two with the *origin 380 381 try { 381 Orthovector1 =Plane(TopOrigin-> getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();382 Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal(); 382 383 } 383 384 catch(LinearDependenceException &excp){ … … 400 401 FirstOtherAtom = World::getInstance().createAtom(); 401 402 SecondOtherAtom = World::getInstance().createAtom(); 402 FirstOtherAtom-> setType(1);403 SecondOtherAtom-> setType(1);404 FirstOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity403 FirstOtherAtom->type = elemente->FindElement(1); 404 SecondOtherAtom->type = elemente->FindElement(1); 405 FirstOtherAtom->v = TopReplacement->v; // copy velocity 405 406 FirstOtherAtom->FixedIon = TopReplacement->FixedIon; 406 SecondOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity407 SecondOtherAtom->v = TopReplacement->v; // copy velocity 407 408 SecondOtherAtom->FixedIon = TopReplacement->FixedIon; 408 409 FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father 409 410 SecondOtherAtom->father = NULL; // we are just an added hydrogen with no father 410 bondangle = TopOrigin-> getType()->HBondAngle[1];411 bondangle = TopOrigin->type->HBondAngle[1]; 411 412 if (bondangle == -1) { 412 413 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl); … … 422 423 // Log() << Verbose(0) << endl; 423 424 // Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl; 424 FirstOtherAtom-> Zero();425 SecondOtherAtom-> Zero();425 FirstOtherAtom->x.Zero(); 426 SecondOtherAtom->x.Zero(); 426 427 for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction) 427 FirstOtherAtom-> set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));428 SecondOtherAtom-> set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));429 } 430 FirstOtherAtom-> Scale(BondRescale); // rescale by correct BondDistance431 SecondOtherAtom-> Scale(BondRescale);428 FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)); 429 SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)); 430 } 431 FirstOtherAtom->x *= BondRescale; // rescale by correct BondDistance 432 SecondOtherAtom->x *= BondRescale; 432 433 //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl; 433 *FirstOtherAtom += TopOrigin->getPosition(); 434 *SecondOtherAtom += TopOrigin->getPosition(); 434 for(int i=NDIM;i--;) { // and make relative to origin atom 435 FirstOtherAtom->x[i] += TopOrigin->x[i]; 436 SecondOtherAtom->x[i] += TopOrigin->x[i]; 437 } 435 438 // ... and add to molecule 436 439 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom); … … 454 457 SecondOtherAtom = World::getInstance().createAtom(); 455 458 ThirdOtherAtom = World::getInstance().createAtom(); 456 FirstOtherAtom-> setType(1);457 SecondOtherAtom-> setType(1);458 ThirdOtherAtom-> setType(1);459 FirstOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity459 FirstOtherAtom->type = elemente->FindElement(1); 460 SecondOtherAtom->type = elemente->FindElement(1); 461 ThirdOtherAtom->type = elemente->FindElement(1); 462 FirstOtherAtom->v = TopReplacement->v; // copy velocity 460 463 FirstOtherAtom->FixedIon = TopReplacement->FixedIon; 461 SecondOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity464 SecondOtherAtom->v = TopReplacement->v; // copy velocity 462 465 SecondOtherAtom->FixedIon = TopReplacement->FixedIon; 463 ThirdOtherAtom-> AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity466 ThirdOtherAtom->v = TopReplacement->v; // copy velocity 464 467 ThirdOtherAtom->FixedIon = TopReplacement->FixedIon; 465 468 FirstOtherAtom->father = NULL; // we are just an added hydrogen with no father … … 484 487 485 488 // create correct coordination for the three atoms 486 alpha = (TopOrigin-> getType()->HBondAngle[2])/180.*M_PI/2.; // retrieve triple bond angle from database489 alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.; // retrieve triple bond angle from database 487 490 l = BondRescale; // desired bond length 488 491 b = 2.*l*sin(alpha); // base length of isosceles triangle … … 495 498 factors[1] = f; 496 499 factors[2] = 0.; 497 FirstOtherAtom-> LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);500 FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors); 498 501 factors[1] = -0.5*f; 499 502 factors[2] = g; 500 SecondOtherAtom-> LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);503 SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors); 501 504 factors[2] = -g; 502 ThirdOtherAtom-> LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);505 ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors); 503 506 504 507 // rescale each to correct BondDistance … … 508 511 509 512 // and relative to *origin atom 510 *FirstOtherAtom += TopOrigin->getPosition();511 *SecondOtherAtom += TopOrigin->getPosition();512 *ThirdOtherAtom += TopOrigin->getPosition();513 FirstOtherAtom->x += TopOrigin->x; 514 SecondOtherAtom->x += TopOrigin->x; 515 ThirdOtherAtom->x += TopOrigin->x; 513 516 514 517 // ... and add to molecule … … 586 589 *item >> x[1]; 587 590 *item >> x[2]; 588 Walker-> setType(elemente->FindElement(shorthand));589 if (Walker-> getType()== NULL) {591 Walker->type = elemente->FindElement(shorthand); 592 if (Walker->type == NULL) { 590 593 DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H."); 591 Walker-> setType(1);594 Walker->type = elemente->FindElement(1); 592 595 } 593 596 if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) { … … 596 599 Walker->Trajectory.F.resize(MDSteps+10); 597 600 } 598 Walker->setPosition(Vector(x));599 601 for(j=NDIM;j--;) { 602 Walker->x[j] = x[j]; 600 603 Walker->Trajectory.R.at(MDSteps-1)[j] = x[j]; 601 604 Walker->Trajectory.U.at(MDSteps-1)[j] = 0; … … 691 694 atom1->RegisterBond(Binder); 692 695 atom2->RegisterBond(Binder); 693 if ((atom1-> getType() != NULL) && (atom1->getType()->Z != 1) && (atom2->getType() != NULL) && (atom2->getType()->Z != 1))696 if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1)) 694 697 NoNonBonds++; 695 698 … … 765 768 ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom()."); 766 769 OBSERVE; 767 formula-=pointer-> getType();770 formula-=pointer->type; 768 771 RemoveBonds(pointer); 769 772 erase(pointer); … … 779 782 if (pointer == NULL) 780 783 return false; 781 formula-=pointer-> getType();784 formula-=pointer->type; 782 785 erase(pointer); 783 786 return true; … … 966 969 for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) { 967 970 (*iter)->nr = i; // update number in molecule (for easier referencing in FragmentMolecule lateron) 968 if ((*iter)-> getType()->Z != 1) // count non-hydrogen atoms whilst at it971 if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it 969 972 NoNonHydrogen++; 970 973 stringstream sstr; 971 sstr << (*iter)-> getType()->symbol << (*iter)->nr+1;974 sstr << (*iter)->type->symbol << (*iter)->nr+1; 972 975 (*iter)->setName(sstr.str()); 973 976 DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
Note:
See TracChangeset
for help on using the changeset viewer.
