Changeset 833b15 for src/molecule_geometry.cpp
- Timestamp:
- Sep 9, 2014, 7:42:32 PM (10 years ago)
- Branches:
- 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
- Children:
- a090e3
- Parents:
- 2d701e
- git-author:
- Frederik Heber <heber@…> (09/01/14 15:54:02)
- git-committer:
- Frederik Heber <heber@…> (09/09/14 19:42:32)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/molecule_geometry.cpp
r2d701e r833b15 58 58 /************************************* Functions for class molecule *********************************/ 59 59 60 /** Returns vector pointing to center of the domain. 61 * \return pointer to center of the domain 62 */ 63 #ifdef HAVE_INLINE 64 inline 65 #else 66 static 67 #endif 68 const Vector DetermineCenterOfBox() 69 { 70 Vector a(0.5,0.5,0.5); 71 const RealSpaceMatrix &M = World::getInstance().getDomain().getM(); 72 a *= M; 73 return a; 74 } 60 75 61 76 /** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths. … … 65 80 { 66 81 bool status = true; 67 const Vector *Center = DetermineCenterOfAll();68 const Vector *CenterBox = DetermineCenterOfBox();82 const Vector Center = DetermineCenterOfAll(); 83 const Vector CenterBox = DetermineCenterOfBox(); 69 84 Box &domain = World::getInstance().getDomain(); 70 85 71 86 // go through all atoms 72 const Vector difference = *CenterBox - *Center; 73 Translate(&difference); 87 Translate(CenterBox - Center); 74 88 getAtomSet().transformNodes(boost::bind(&Box::enforceBoundaryConditions,domain,_1)); 75 89 76 delete(Center);77 delete(CenterBox);78 90 return status; 79 } ;91 } 80 92 81 93 … … 92 104 93 105 return status; 94 } ;106 } 95 107 96 108 /** Centers the edge of the atoms at (0,0,0). 97 * \param *out output stream for debugging 98 * \param *max coordinates of other edge, specifying box dimensions. 99 */ 100 void molecule::CenterEdge(Vector *max) 109 */ 110 void molecule::CenterEdge() 101 111 { 102 112 const_iterator iter = begin(); 103 113 if (iter != end()) { //list not empty? 104 114 Vector min = (*begin())->getPosition(); 105 *max = min;106 115 for (;iter != end(); ++iter) { // continue with second if present 107 116 const Vector ¤tPos = (*iter)->getPosition(); 108 for (size_t i=0;i<NDIM;++i) {117 for (size_t i=0;i<NDIM;++i) 109 118 if (min[i] > currentPos[i]) 110 119 min[i] = currentPos[i]; 111 if ((*max)[i] < currentPos[i]) 112 (*max)[i] = currentPos[i]; 113 } 114 } 115 LOG(1, "INFO: Maximum is " << *max << ", Minimum is " << min << "."); 116 const Vector temp = -1.*min; 117 Translate(&temp); 118 } 119 }; 120 } 121 Translate(-1.*min); 122 } 123 } 120 124 121 125 /** Centers the center of the atoms at (0,0,0). … … 136 140 } 137 141 Center.Scale(-1./(double)Num); // divide through total number (and sign for direction) 138 Translate( &Center);139 } 140 } ;142 Translate(Center); 143 } 144 } 141 145 142 146 /** Returns vector pointing to center of all atoms. 143 147 * \return pointer to center of all vector 144 148 */ 145 Vector *molecule::DetermineCenterOfAll() const149 const Vector molecule::DetermineCenterOfAll() const 146 150 { 147 151 const_iterator iter = begin(); // start at first in list 148 Vector *a = new Vector();152 Vector a; 149 153 double Num = 0; 150 154 151 a ->Zero();155 a.Zero(); 152 156 153 157 if (iter != end()) { //list not empty? 154 158 for (; iter != end(); ++iter) { // continue with second if present 155 159 Num++; 156 (*a)+= (*iter)->getPosition();157 } 158 a ->Scale(1./(double)Num); // divide through total mass (and sign for direction)160 a += (*iter)->getPosition(); 161 } 162 a.Scale(1./(double)Num); // divide through total mass (and sign for direction) 159 163 } 160 164 return a; 161 }; 162 163 /** Returns vector pointing to center of the domain. 164 * \return pointer to center of the domain 165 */ 166 Vector * molecule::DetermineCenterOfBox() const 167 { 168 Vector *a = new Vector(0.5,0.5,0.5); 169 const RealSpaceMatrix &M = World::getInstance().getDomain().getM(); 170 (*a) *= M; 171 return a; 172 }; 165 } 166 173 167 174 168 /** Returns vector pointing to center of gravity. … … 176 170 * \return pointer to center of gravity vector 177 171 */ 178 Vector *molecule::DetermineCenterOfGravity() const172 const Vector molecule::DetermineCenterOfGravity() const 179 173 { 180 174 const_iterator iter = begin(); // start at first in list 181 Vector *a = new Vector();175 Vector a; 182 176 Vector tmp; 183 177 double Num = 0; 184 178 185 a ->Zero();179 a.Zero(); 186 180 187 181 if (iter != end()) { //list not empty? … … 189 183 Num += (*iter)->getType()->getMass(); 190 184 tmp = (*iter)->getType()->getMass() * (*iter)->getPosition(); 191 (*a)+= tmp;192 } 193 a ->Scale(1./Num); // divide through total mass194 } 195 LOG(1, "INFO: Resulting center of gravity: " << *a << ".");185 a += tmp; 186 } 187 a.Scale(1./Num); // divide through total mass 188 } 189 LOG(1, "INFO: Resulting center of gravity: " << a << "."); 196 190 return a; 197 } ;191 } 198 192 199 193 /** Centers the center of gravity of the atoms at (0,0,0). … … 205 199 Vector NewCenter; 206 200 DeterminePeriodicCenter(NewCenter); 207 NewCenter *= -1.; 208 Translate(&NewCenter); 209 }; 201 Translate(-1.*NewCenter); 202 } 210 203 211 204 … … 214 207 * \param *center return vector for translation vector 215 208 */ 216 void molecule::CenterAtVector(Vector *newcenter) 217 { 218 const Vector temp = -1.**newcenter; 219 Translate(&temp); 220 }; 209 void molecule::CenterAtVector(const Vector &newcenter) 210 { 211 Translate(-1.*newcenter); 212 } 221 213 222 214 /** Calculate the inertia tensor of a the molecule. … … 227 219 { 228 220 RealSpaceMatrix InertiaTensor; 229 Vector *CenterOfGravity = DetermineCenterOfGravity();221 const Vector CenterOfGravity = DetermineCenterOfGravity(); 230 222 231 223 // reset inertia tensor … … 235 227 for (const_iterator iter = begin(); iter != end(); ++iter) { 236 228 Vector x = (*iter)->getPosition(); 237 x -= *CenterOfGravity;229 x -= CenterOfGravity; 238 230 const double mass = (*iter)->getType()->getMass(); 239 231 InertiaTensor.at(0,0) += mass*(x[1]*x[1] + x[2]*x[2]); … … 250 242 LOG(1, "INFO: The inertia tensor of molecule " << getName() << " is:" << InertiaTensor); 251 243 252 delete CenterOfGravity;253 244 return InertiaTensor; 254 245 } … … 261 252 void molecule::RotateToPrincipalAxisSystem(const Vector &Axis) 262 253 { 263 Vector *CenterOfGravity = DetermineCenterOfGravity();254 const Vector CenterOfGravity = DetermineCenterOfGravity(); 264 255 RealSpaceMatrix InertiaTensor = getInertiaTensor(); 265 256 … … 273 264 274 265 // obtain first column, eigenvector to biggest eigenvalue 275 Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent()));266 const Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent())); 276 267 Vector DesiredAxis(Axis.getNormalized()); 277 268 … … 287 278 // and rotate 288 279 for (iterator iter = begin(); iter != end(); ++iter) { 289 *(*iter) -= *CenterOfGravity;280 *(*iter) -= CenterOfGravity; 290 281 (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), alpha)); 291 *(*iter) += *CenterOfGravity;282 *(*iter) += CenterOfGravity; 292 283 } 293 284 LOG(0, "STATUS: done."); 294 295 delete CenterOfGravity;296 285 } 297 286 … … 304 293 * x=(**factor) * x (as suggested by comment) 305 294 */ 306 void molecule::Scale(const double * * constfactor)295 void molecule::Scale(const double *factor) 307 296 { 308 297 for (iterator iter = begin(); iter != end(); ++iter) { 309 298 for (size_t j=0;j<(*iter)->getTrajectorySize();j++) { 310 299 Vector temp = (*iter)->getPositionAtStep(j); 311 temp.ScaleAll( *factor);300 temp.ScaleAll(factor); 312 301 (*iter)->setPositionAtStep(j,temp); 313 302 } … … 318 307 * \param trans[] translation vector. 319 308 */ 320 void molecule::Translate(const Vector *trans)321 { 322 getAtomSet().translate( *trans);309 void molecule::Translate(const Vector &trans) 310 { 311 getAtomSet().translate(trans); 323 312 }; 324 313 … … 327 316 * TODO treatment of trajectories missing 328 317 */ 329 void molecule::TranslatePeriodically(const Vector *trans)318 void molecule::TranslatePeriodically(const Vector &trans) 330 319 { 331 320 Translate(trans); … … 338 327 * \param n[] normal vector of mirror plane. 339 328 */ 340 void molecule::Mirror(const Vector *n)341 { 342 Plane p( *n,0);329 void molecule::Mirror(const Vector &n) 330 { 331 Plane p(n,0); 343 332 getAtomSet().transformNodes(boost::bind(&Plane::mirrorVector,p,_1)); 344 333 }; … … 356 345 Vector Testvector, Translationvector; 357 346 Vector Center; 358 BondGraph *BG = World::getInstance().getBondGraph();347 const BondGraph * const BG = World::getInstance().getBondGraph(); 359 348 360 349 do { … … 407 396 408 397 Center.Scale(1./static_cast<double>(getAtomCount())); 409 CenterAtVector( &Center);398 CenterAtVector(Center); 410 399 }; 411 400 … … 413 402 * \param n[] alignment vector. 414 403 */ 415 void molecule::Align( Vector *n)404 void molecule::Align(const Vector &n) 416 405 { 417 406 double alpha, tmp; 418 407 Vector z_axis; 408 Vector alignment(n); 419 409 z_axis[0] = 0.; 420 410 z_axis[1] = 0.; … … 423 413 // rotate on z-x plane 424 414 LOG(0, "Begin of Aligning all atoms."); 425 alpha = atan(- n->at(0)/n->at(2));415 alpha = atan(-alignment.at(0)/alignment.at(2)); 426 416 LOG(1, "INFO: Z-X-angle: " << alpha << " ... "); 427 417 for (iterator iter = begin(); iter != end(); ++iter) { … … 437 427 } 438 428 // rotate n vector 439 tmp = n->at(0);440 n->at(0) = cos(alpha) * tmp + sin(alpha) * n->at(2);441 n->at(2) = -sin(alpha) * tmp + cos(alpha) * n->at(2);442 LOG(1, "alignment vector after first rotation: " << n);429 tmp = alignment.at(0); 430 alignment.at(0) = cos(alpha) * tmp + sin(alpha) * alignment.at(2); 431 alignment.at(2) = -sin(alpha) * tmp + cos(alpha) * alignment.at(2); 432 LOG(1, "alignment vector after first rotation: " << alignment); 443 433 444 434 // rotate on z-y plane 445 alpha = atan(- n->at(1)/n->at(2));435 alpha = atan(-alignment.at(1)/alignment.at(2)); 446 436 LOG(1, "INFO: Z-Y-angle: " << alpha << " ... "); 447 437 for (iterator iter = begin(); iter != end(); ++iter) { … … 457 447 } 458 448 // rotate n vector (for consistency check) 459 tmp = n->at(1); 460 n->at(1) = cos(alpha) * tmp + sin(alpha) * n->at(2); 461 n->at(2) = -sin(alpha) * tmp + cos(alpha) * n->at(2); 462 463 464 LOG(1, "alignment vector after second rotation: " << n); 449 tmp = alignment.at(1); 450 alignment.at(1) = cos(alpha) * tmp + sin(alpha) * alignment.at(2); 451 alignment.at(2) = -sin(alpha) * tmp + cos(alpha) * alignment.at(2); 452 453 LOG(1, "alignment vector after second rotation: " << alignment); 465 454 LOG(0, "End of Aligning all atoms."); 466 455 };
Note:
See TracChangeset
for help on using the changeset viewer.