source: src/vector.cpp@ 8c54a3

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 8c54a3 was 9c20aa, checked in by Saskia Metzler <metzler@…>, 15 years ago

new function isNull(), output operator takes const argument

  • Property mode set to 100755
File size: 25.5 KB
Line 
1/** \file vector.cpp
2 *
3 * Function implementations for the class vector.
4 *
5 */
6
7#include "molecules.hpp"
8
9
10/************************************ Functions for class vector ************************************/
11
12/** Constructor of class vector.
13 */
14Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
15
16/** Constructor of class vector.
17 */
18Vector::Vector(double x1, double x2, double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
19
20/** Desctructor of class vector.
21 */
22Vector::~Vector() {};
23
24/** Calculates square of distance between this and another vector.
25 * \param *y array to second vector
26 * \return \f$| x - y |^2\f$
27 */
28double Vector::DistanceSquared(const Vector *y) const
29{
30 double res = 0.;
31 for (int i=NDIM;i--;)
32 res += (x[i]-y->x[i])*(x[i]-y->x[i]);
33 return (res);
34};
35
36/** Calculates distance between this and another vector.
37 * \param *y array to second vector
38 * \return \f$| x - y |\f$
39 */
40double Vector::Distance(const Vector *y) const
41{
42 double res = 0.;
43 for (int i=NDIM;i--;)
44 res += (x[i]-y->x[i])*(x[i]-y->x[i]);
45 return (sqrt(res));
46};
47
48/** Calculates distance between this and another vector in a periodic cell.
49 * \param *y array to second vector
50 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
51 * \return \f$| x - y |\f$
52 */
53double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
54{
55 double res = Distance(y), tmp, matrix[NDIM*NDIM];
56 Vector Shiftedy, TranslationVector;
57 int N[NDIM];
58 matrix[0] = cell_size[0];
59 matrix[1] = cell_size[1];
60 matrix[2] = cell_size[3];
61 matrix[3] = cell_size[1];
62 matrix[4] = cell_size[2];
63 matrix[5] = cell_size[4];
64 matrix[6] = cell_size[3];
65 matrix[7] = cell_size[4];
66 matrix[8] = cell_size[5];
67 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
68 for (N[0]=-1;N[0]<=1;N[0]++)
69 for (N[1]=-1;N[1]<=1;N[1]++)
70 for (N[2]=-1;N[2]<=1;N[2]++) {
71 // create the translation vector
72 TranslationVector.Zero();
73 for (int i=NDIM;i--;)
74 TranslationVector.x[i] = (double)N[i];
75 TranslationVector.MatrixMultiplication(matrix);
76 // add onto the original vector to compare with
77 Shiftedy.CopyVector(y);
78 Shiftedy.AddVector(&TranslationVector);
79 // get distance and compare with minimum so far
80 tmp = Distance(&Shiftedy);
81 if (tmp < res) res = tmp;
82 }
83 return (res);
84};
85
86/** Calculates distance between this and another vector in a periodic cell.
87 * \param *y array to second vector
88 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
89 * \return \f$| x - y |^2\f$
90 */
91double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
92{
93 double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
94 Vector Shiftedy, TranslationVector;
95 int N[NDIM];
96 matrix[0] = cell_size[0];
97 matrix[1] = cell_size[1];
98 matrix[2] = cell_size[3];
99 matrix[3] = cell_size[1];
100 matrix[4] = cell_size[2];
101 matrix[5] = cell_size[4];
102 matrix[6] = cell_size[3];
103 matrix[7] = cell_size[4];
104 matrix[8] = cell_size[5];
105 // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
106 for (N[0]=-1;N[0]<=1;N[0]++)
107 for (N[1]=-1;N[1]<=1;N[1]++)
108 for (N[2]=-1;N[2]<=1;N[2]++) {
109 // create the translation vector
110 TranslationVector.Zero();
111 for (int i=NDIM;i--;)
112 TranslationVector.x[i] = (double)N[i];
113 TranslationVector.MatrixMultiplication(matrix);
114 // add onto the original vector to compare with
115 Shiftedy.CopyVector(y);
116 Shiftedy.AddVector(&TranslationVector);
117 // get distance and compare with minimum so far
118 tmp = DistanceSquared(&Shiftedy);
119 if (tmp < res) res = tmp;
120 }
121 return (res);
122};
123
124/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
125 * \param *out ofstream for debugging messages
126 * Tries to translate a vector into each adjacent neighbouring cell.
127 */
128void Vector::KeepPeriodic(ofstream *out, double *matrix)
129{
130// int N[NDIM];
131// bool flag = false;
132 //vector Shifted, TranslationVector;
133 Vector TestVector;
134// *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
135// *out << Verbose(2) << "Vector is: ";
136// Output(out);
137// *out << endl;
138 TestVector.CopyVector(this);
139 TestVector.InverseMatrixMultiplication(matrix);
140 for(int i=NDIM;i--;) { // correct periodically
141 if (TestVector.x[i] < 0) { // get every coefficient into the interval [0,1)
142 TestVector.x[i] += ceil(TestVector.x[i]);
143 } else {
144 TestVector.x[i] -= floor(TestVector.x[i]);
145 }
146 }
147 TestVector.MatrixMultiplication(matrix);
148 CopyVector(&TestVector);
149// *out << Verbose(2) << "New corrected vector is: ";
150// Output(out);
151// *out << endl;
152// *out << Verbose(1) << "End of KeepPeriodic." << endl;
153};
154
155/** Calculates scalar product between this and another vector.
156 * \param *y array to second vector
157 * \return \f$\langle x, y \rangle\f$
158 */
159double Vector::ScalarProduct(const Vector *y) const
160{
161 double res = 0.;
162 for (int i=NDIM;i--;)
163 res += x[i]*y->x[i];
164 return (res);
165};
166
167
168/** Calculates VectorProduct between this and another vector.
169 * -# returns the Product in place of vector from which it was initiated
170 * -# ATTENTION: Only three dim.
171 * \param *y array to vector with which to calculate crossproduct
172 * \return \f$ x \times y \f&
173 */
174void Vector::VectorProduct(const Vector *y)
175{
176 Vector tmp;
177 tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
178 tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
179 tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
180 this->CopyVector(&tmp);
181
182};
183
184
185/** projects this vector onto plane defined by \a *y.
186 * \param *y normal vector of plane
187 * \return \f$\langle x, y \rangle\f$
188 */
189void Vector::ProjectOntoPlane(const Vector *y)
190{
191 Vector tmp;
192 tmp.CopyVector(y);
193 tmp.Normalize();
194 tmp.Scale(ScalarProduct(&tmp));
195 this->SubtractVector(&tmp);
196};
197
198/** Calculates the projection of a vector onto another \a *y.
199 * \param *y array to second vector
200 * \return \f$\langle x, y \rangle\f$
201 */
202double Vector::Projection(const Vector *y) const
203{
204 return (ScalarProduct(y));
205};
206
207/** Calculates norm of this vector.
208 * \return \f$|x|\f$
209 */
210double Vector::Norm() const
211{
212 double res = 0.;
213 for (int i=NDIM;i--;)
214 res += this->x[i]*this->x[i];
215 return (sqrt(res));
216};
217
218/** Normalizes this vector.
219 */
220void Vector::Normalize()
221{
222 double res = 0.;
223 for (int i=NDIM;i--;)
224 res += this->x[i]*this->x[i];
225 if (fabs(res) > MYEPSILON)
226 res = 1./sqrt(res);
227 Scale(&res);
228};
229
230/** Zeros all components of this vector.
231 */
232void Vector::Zero()
233{
234 for (int i=NDIM;i--;)
235 this->x[i] = 0.;
236};
237
238/** Zeros all components of this vector.
239 */
240void Vector::One(double one)
241{
242 for (int i=NDIM;i--;)
243 this->x[i] = one;
244};
245
246/** Initialises all components of this vector.
247 */
248void Vector::Init(double x1, double x2, double x3)
249{
250 x[0] = x1;
251 x[1] = x2;
252 x[2] = x3;
253};
254
255/** Checks whether vector has all components zero.
256 * @return true - vector is zero, false - vector is not
257 */
258bool Vector::IsNull()
259{
260 return (fabs(x[0]+x[1]+x[2]) < MYEPSILON);
261};
262
263/** Calculates the angle between this and another vector.
264 * \param *y array to second vector
265 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
266 */
267double Vector::Angle(const Vector *y) const
268{
269 double angle = this->ScalarProduct(y)/Norm()/y->Norm();
270 // -1-MYEPSILON occured due to numerical imprecision, catch ...
271 //cout << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
272 if (angle < -1)
273 angle = -1;
274 if (angle > 1)
275 angle = 1;
276 return acos(angle);
277};
278
279/** Rotates the vector around the axis given by \a *axis by an angle of \a alpha.
280 * \param *axis rotation axis
281 * \param alpha rotation angle in radian
282 */
283void Vector::RotateVector(const Vector *axis, const double alpha)
284{
285 Vector a,y;
286 // normalise this vector with respect to axis
287 a.CopyVector(this);
288 a.Scale(Projection(axis));
289 SubtractVector(&a);
290 // construct normal vector
291 y.MakeNormalVector(axis,this);
292 y.Scale(Norm());
293 // scale normal vector by sine and this vector by cosine
294 y.Scale(sin(alpha));
295 Scale(cos(alpha));
296 // add scaled normal vector onto this vector
297 AddVector(&y);
298 // add part in axis direction
299 AddVector(&a);
300};
301
302/** Sums vector \a to this lhs component-wise.
303 * \param a base vector
304 * \param b vector components to add
305 * \return lhs + a
306 */
307Vector& operator+=(Vector& a, const Vector& b)
308{
309 a.AddVector(&b);
310 return a;
311};
312/** factor each component of \a a times a double \a m.
313 * \param a base vector
314 * \param m factor
315 * \return lhs.x[i] * m
316 */
317Vector& operator*=(Vector& a, const double m)
318{
319 a.Scale(m);
320 return a;
321};
322
323/** Sums two vectors \a and \b component-wise.
324 * \param a first vector
325 * \param b second vector
326 * \return a + b
327 */
328Vector& operator+(const Vector& a, const Vector& b)
329{
330 Vector *x = new Vector;
331 x->CopyVector(&a);
332 x->AddVector(&b);
333 return *x;
334};
335
336/** Factors given vector \a a times \a m.
337 * \param a vector
338 * \param m factor
339 * \return a + b
340 */
341Vector& operator*(const Vector& a, const double m)
342{
343 Vector *x = new Vector;
344 x->CopyVector(&a);
345 x->Scale(m);
346 return *x;
347};
348
349/** Prints a 3dim vector.
350 * prints no end of line.
351 * \param *out output stream
352 */
353bool Vector::Output(ofstream *out) const
354{
355 if (out != NULL) {
356 *out << "(";
357 for (int i=0;i<NDIM;i++) {
358 *out << x[i];
359 if (i != 2)
360 *out << ",";
361 }
362 *out << ")";
363 return true;
364 } else
365 return false;
366};
367
368ostream& operator<<(ostream& ost, const Vector& m)
369{
370 ost << "(";
371 for (int i=0;i<NDIM;i++) {
372 ost << m.x[i];
373 if (i != 2)
374 ost << ",";
375 }
376 ost << ")";
377 return ost;
378};
379
380/** Scales each atom coordinate by an individual \a factor.
381 * \param *factor pointer to scaling factor
382 */
383void Vector::Scale(double **factor)
384{
385 for (int i=NDIM;i--;)
386 x[i] *= (*factor)[i];
387};
388
389void Vector::Scale(double *factor)
390{
391 for (int i=NDIM;i--;)
392 x[i] *= *factor;
393};
394
395void Vector::Scale(double factor)
396{
397 for (int i=NDIM;i--;)
398 x[i] *= factor;
399};
400
401/** Translate atom by given vector.
402 * \param trans[] translation vector.
403 */
404void Vector::Translate(const Vector *trans)
405{
406 for (int i=NDIM;i--;)
407 x[i] += trans->x[i];
408};
409
410/** Do a matrix multiplication.
411 * \param *matrix NDIM_NDIM array
412 */
413void Vector::MatrixMultiplication(double *M)
414{
415 Vector C;
416 // do the matrix multiplication
417 C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
418 C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
419 C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
420 // transfer the result into this
421 for (int i=NDIM;i--;)
422 x[i] = C.x[i];
423};
424
425/** Do a matrix multiplication with \a *matrix' inverse.
426 * \param *matrix NDIM_NDIM array
427 */
428void Vector::InverseMatrixMultiplication(double *A)
429{
430 Vector C;
431 double B[NDIM*NDIM];
432 double detA = RDET3(A);
433 double detAReci;
434
435 // calculate the inverse B
436 if (fabs(detA) > MYEPSILON) {; // RDET3(A) yields precisely zero if A irregular
437 detAReci = 1./detA;
438 B[0] = detAReci*RDET2(A[4],A[5],A[7],A[8]); // A_11
439 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]); // A_12
440 B[2] = detAReci*RDET2(A[1],A[2],A[4],A[5]); // A_13
441 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]); // A_21
442 B[4] = detAReci*RDET2(A[0],A[2],A[6],A[8]); // A_22
443 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]); // A_23
444 B[6] = detAReci*RDET2(A[3],A[4],A[6],A[7]); // A_31
445 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]); // A_32
446 B[8] = detAReci*RDET2(A[0],A[1],A[3],A[4]); // A_33
447
448 // do the matrix multiplication
449 C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
450 C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
451 C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
452 // transfer the result into this
453 for (int i=NDIM;i--;)
454 x[i] = C.x[i];
455 } else {
456 cerr << "ERROR: inverse of matrix does not exists!" << endl;
457 }
458};
459
460
461/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
462 * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
463 * \param *x1 first vector
464 * \param *x2 second vector
465 * \param *x3 third vector
466 * \param *factors three-component vector with the factor for each given vector
467 */
468void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
469{
470 for(int i=NDIM;i--;)
471 x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
472};
473
474/** Mirrors atom against a given plane.
475 * \param n[] normal vector of mirror plane.
476 */
477void Vector::Mirror(const Vector *n)
478{
479 double projection;
480 projection = ScalarProduct(n)/n->ScalarProduct(n); // remove constancy from n (keep as logical one)
481 // withdraw projected vector twice from original one
482 cout << Verbose(1) << "Vector: ";
483 Output((ofstream *)&cout);
484 cout << "\t";
485 for (int i=NDIM;i--;)
486 x[i] -= 2.*projection*n->x[i];
487 cout << "Projected vector: ";
488 Output((ofstream *)&cout);
489 cout << endl;
490};
491
492/** Calculates normal vector for three given vectors (being three points in space).
493 * Makes this vector orthonormal to the three given points, making up a place in 3d space.
494 * \param *y1 first vector
495 * \param *y2 second vector
496 * \param *y3 third vector
497 * \return true - success, vectors are linear independent, false - failure due to linear dependency
498 */
499bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
500{
501 Vector x1, x2;
502
503 x1.CopyVector(y1);
504 x1.SubtractVector(y2);
505 x2.CopyVector(y3);
506 x2.SubtractVector(y2);
507 if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
508 cout << Verbose(4) << "Given vectors are linear dependent." << endl;
509 return false;
510 }
511// cout << Verbose(4) << "relative, first plane coordinates:";
512// x1.Output((ofstream *)&cout);
513// cout << endl;
514// cout << Verbose(4) << "second plane coordinates:";
515// x2.Output((ofstream *)&cout);
516// cout << endl;
517
518 this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
519 this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
520 this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
521 Normalize();
522
523 return true;
524};
525
526
527/** Calculates orthonormal vector to two given vectors.
528 * Makes this vector orthonormal to two given vectors. This is very similar to the other
529 * vector::MakeNormalVector(), only there three points whereas here two difference
530 * vectors are given.
531 * \param *x1 first vector
532 * \param *x2 second vector
533 * \return true - success, vectors are linear independent, false - failure due to linear dependency
534 */
535bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
536{
537 Vector x1,x2;
538 x1.CopyVector(y1);
539 x2.CopyVector(y2);
540 Zero();
541 if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
542 cout << Verbose(4) << "Given vectors are linear dependent." << endl;
543 return false;
544 }
545// cout << Verbose(4) << "relative, first plane coordinates:";
546// x1.Output((ofstream *)&cout);
547// cout << endl;
548// cout << Verbose(4) << "second plane coordinates:";
549// x2.Output((ofstream *)&cout);
550// cout << endl;
551
552 this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
553 this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
554 this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
555 Normalize();
556
557 return true;
558};
559
560/** Calculates orthonormal vector to one given vectors.
561 * Just subtracts the projection onto the given vector from this vector.
562 * \param *x1 vector
563 * \return true - success, false - vector is zero
564 */
565bool Vector::MakeNormalVector(const Vector *y1)
566{
567 bool result = false;
568 Vector x1;
569 x1.CopyVector(y1);
570 x1.Scale(x1.Projection(this));
571 SubtractVector(&x1);
572 for (int i=NDIM;i--;)
573 result = result || (fabs(x[i]) > MYEPSILON);
574
575 return result;
576};
577
578/** Creates this vector as one of the possible orthonormal ones to the given one.
579 * Just scan how many components of given *vector are unequal to zero and
580 * try to get the skp of both to be zero accordingly.
581 * \param *vector given vector
582 * \return true - success, false - failure (null vector given)
583 */
584bool Vector::GetOneNormalVector(const Vector *GivenVector)
585{
586 int Components[NDIM]; // contains indices of non-zero components
587 int Last = 0; // count the number of non-zero entries in vector
588 int j; // loop variables
589 double norm;
590
591 cout << Verbose(4);
592 GivenVector->Output((ofstream *)&cout);
593 cout << endl;
594 for (j=NDIM;j--;)
595 Components[j] = -1;
596 // find two components != 0
597 for (j=0;j<NDIM;j++)
598 if (fabs(GivenVector->x[j]) > MYEPSILON)
599 Components[Last++] = j;
600 cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
601
602 switch(Last) {
603 case 3: // threecomponent system
604 case 2: // two component system
605 norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
606 x[Components[2]] = 0.;
607 // in skp both remaining parts shall become zero but with opposite sign and third is zero
608 x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
609 x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
610 return true;
611 break;
612 case 1: // one component system
613 // set sole non-zero component to 0, and one of the other zero component pendants to 1
614 x[(Components[0]+2)%NDIM] = 0.;
615 x[(Components[0]+1)%NDIM] = 1.;
616 x[Components[0]] = 0.;
617 return true;
618 break;
619 default:
620 return false;
621 }
622};
623
624/** Determines paramter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
625 * \param *A first plane vector
626 * \param *B second plane vector
627 * \param *C third plane vector
628 * \return scaling parameter for this vector
629 */
630double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
631{
632// cout << Verbose(3) << "For comparison: ";
633// cout << "A " << A->Projection(this) << "\t";
634// cout << "B " << B->Projection(this) << "\t";
635// cout << "C " << C->Projection(this) << "\t";
636// cout << endl;
637 return A->Projection(this);
638};
639
640/** Creates a new vector as the one with least square distance to a given set of \a vectors.
641 * \param *vectors set of vectors
642 * \param num number of vectors
643 * \return true if success, false if failed due to linear dependency
644 */
645bool Vector::LSQdistance(Vector **vectors, int num)
646{
647 int j;
648
649 for (j=0;j<num;j++) {
650 cout << Verbose(1) << j << "th atom's vector: ";
651 (vectors[j])->Output((ofstream *)&cout);
652 cout << endl;
653 }
654
655 int np = 3;
656 struct LSQ_params par;
657
658 const gsl_multimin_fminimizer_type *T =
659 gsl_multimin_fminimizer_nmsimplex;
660 gsl_multimin_fminimizer *s = NULL;
661 gsl_vector *ss, *y;
662 gsl_multimin_function minex_func;
663
664 size_t iter = 0, i;
665 int status;
666 double size;
667
668 /* Initial vertex size vector */
669 ss = gsl_vector_alloc (np);
670 y = gsl_vector_alloc (np);
671
672 /* Set all step sizes to 1 */
673 gsl_vector_set_all (ss, 1.0);
674
675 /* Starting point */
676 par.vectors = vectors;
677 par.num = num;
678
679 for (i=NDIM;i--;)
680 gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
681
682 /* Initialize method and iterate */
683 minex_func.f = &LSQ;
684 minex_func.n = np;
685 minex_func.params = (void *)&par;
686
687 s = gsl_multimin_fminimizer_alloc (T, np);
688 gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
689
690 do
691 {
692 iter++;
693 status = gsl_multimin_fminimizer_iterate(s);
694
695 if (status)
696 break;
697
698 size = gsl_multimin_fminimizer_size (s);
699 status = gsl_multimin_test_size (size, 1e-2);
700
701 if (status == GSL_SUCCESS)
702 {
703 printf ("converged to minimum at\n");
704 }
705
706 printf ("%5d ", (int)iter);
707 for (i = 0; i < (size_t)np; i++)
708 {
709 printf ("%10.3e ", gsl_vector_get (s->x, i));
710 }
711 printf ("f() = %7.3f size = %.3f\n", s->fval, size);
712 }
713 while (status == GSL_CONTINUE && iter < 100);
714
715 for (i=(size_t)np;i--;)
716 this->x[i] = gsl_vector_get(s->x, i);
717 gsl_vector_free(y);
718 gsl_vector_free(ss);
719 gsl_multimin_fminimizer_free (s);
720
721 return true;
722};
723
724/** Adds vector \a *y componentwise.
725 * \param *y vector
726 */
727void Vector::AddVector(const Vector *y)
728{
729 for (int i=NDIM;i--;)
730 this->x[i] += y->x[i];
731}
732
733/** Adds vector \a *y componentwise.
734 * \param *y vector
735 */
736void Vector::SubtractVector(const Vector *y)
737{
738 for (int i=NDIM;i--;)
739 this->x[i] -= y->x[i];
740}
741
742/** Copy vector \a *y componentwise.
743 * \param *y vector
744 */
745void Vector::CopyVector(const Vector *y)
746{
747 for (int i=NDIM;i--;)
748 this->x[i] = y->x[i];
749}
750
751
752/** Asks for position, checks for boundary.
753 * \param cell_size unitary size of cubic cell, coordinates must be within 0...cell_size
754 * \param check whether bounds shall be checked (true) or not (false)
755 */
756void Vector::AskPosition(double *cell_size, bool check)
757{
758 char coords[3] = {'x','y','z'};
759 int j = -1;
760 for (int i=0;i<3;i++) {
761 j += i+1;
762 do {
763 cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
764 cin >> x[i];
765 } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
766 }
767};
768
769/** Solves a vectorial system consisting of two orthogonal statements and a norm statement.
770 * This is linear system of equations to be solved, however of the three given (skp of this vector\
771 * with either of the three hast to be zero) only two are linear independent. The third equation
772 * is that the vector should be of magnitude 1 (orthonormal). This all leads to a case-based solution
773 * where very often it has to be checked whether a certain value is zero or not and thus forked into
774 * another case.
775 * \param *x1 first vector
776 * \param *x2 second vector
777 * \param *y third vector
778 * \param alpha first angle
779 * \param beta second angle
780 * \param c norm of final vector
781 * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
782 * \bug this is not yet working properly
783 */
784bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
785{
786 double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
787 double ang; // angle on testing
788 double sign[3];
789 int i,j,k;
790 A = cos(alpha) * x1->Norm() * c;
791 B1 = cos(beta + M_PI/2.) * y->Norm() * c;
792 B2 = cos(beta) * x2->Norm() * c;
793 C = c * c;
794 cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
795 int flag = 0;
796 if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
797 if (fabs(x1->x[1]) > MYEPSILON) {
798 flag = 1;
799 } else if (fabs(x1->x[2]) > MYEPSILON) {
800 flag = 2;
801 } else {
802 return false;
803 }
804 }
805 switch (flag) {
806 default:
807 case 0:
808 break;
809 case 2:
810 flip(&x1->x[0],&x1->x[1]);
811 flip(&x2->x[0],&x2->x[1]);
812 flip(&y->x[0],&y->x[1]);
813 //flip(&x[0],&x[1]);
814 flip(&x1->x[1],&x1->x[2]);
815 flip(&x2->x[1],&x2->x[2]);
816 flip(&y->x[1],&y->x[2]);
817 //flip(&x[1],&x[2]);
818 case 1:
819 flip(&x1->x[0],&x1->x[1]);
820 flip(&x2->x[0],&x2->x[1]);
821 flip(&y->x[0],&y->x[1]);
822 //flip(&x[0],&x[1]);
823 flip(&x1->x[1],&x1->x[2]);
824 flip(&x2->x[1],&x2->x[2]);
825 flip(&y->x[1],&y->x[2]);
826 //flip(&x[1],&x[2]);
827 break;
828 }
829 // now comes the case system
830 D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
831 D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
832 D3 = y->x[0]/x1->x[0]*A-B1;
833 cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
834 if (fabs(D1) < MYEPSILON) {
835 cout << Verbose(2) << "D1 == 0!\n";
836 if (fabs(D2) > MYEPSILON) {
837 cout << Verbose(3) << "D2 != 0!\n";
838 x[2] = -D3/D2;
839 E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
840 E2 = -x1->x[1]/x1->x[0];
841 cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
842 F1 = E1*E1 + 1.;
843 F2 = -E1*E2;
844 F3 = E1*E1 + D3*D3/(D2*D2) - C;
845 cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
846 if (fabs(F1) < MYEPSILON) {
847 cout << Verbose(4) << "F1 == 0!\n";
848 cout << Verbose(4) << "Gleichungssystem linear\n";
849 x[1] = F3/(2.*F2);
850 } else {
851 p = F2/F1;
852 q = p*p - F3/F1;
853 cout << Verbose(4) << "p " << p << "\tq " << q << endl;
854 if (q < 0) {
855 cout << Verbose(4) << "q < 0" << endl;
856 return false;
857 }
858 x[1] = p + sqrt(q);
859 }
860 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
861 } else {
862 cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
863 return false;
864 }
865 } else {
866 E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
867 E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
868 cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
869 F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
870 F2 = -(E1*E2 + D2*D3/(D1*D1));
871 F3 = E1*E1 + D3*D3/(D1*D1) - C;
872 cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
873 if (fabs(F1) < MYEPSILON) {
874 cout << Verbose(3) << "F1 == 0!\n";
875 cout << Verbose(3) << "Gleichungssystem linear\n";
876 x[2] = F3/(2.*F2);
877 } else {
878 p = F2/F1;
879 q = p*p - F3/F1;
880 cout << Verbose(3) << "p " << p << "\tq " << q << endl;
881 if (q < 0) {
882 cout << Verbose(3) << "q < 0" << endl;
883 return false;
884 }
885 x[2] = p + sqrt(q);
886 }
887 x[1] = (-D2 * x[2] - D3)/D1;
888 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
889 }
890 switch (flag) { // back-flipping
891 default:
892 case 0:
893 break;
894 case 2:
895 flip(&x1->x[0],&x1->x[1]);
896 flip(&x2->x[0],&x2->x[1]);
897 flip(&y->x[0],&y->x[1]);
898 flip(&x[0],&x[1]);
899 flip(&x1->x[1],&x1->x[2]);
900 flip(&x2->x[1],&x2->x[2]);
901 flip(&y->x[1],&y->x[2]);
902 flip(&x[1],&x[2]);
903 case 1:
904 flip(&x1->x[0],&x1->x[1]);
905 flip(&x2->x[0],&x2->x[1]);
906 flip(&y->x[0],&y->x[1]);
907 //flip(&x[0],&x[1]);
908 flip(&x1->x[1],&x1->x[2]);
909 flip(&x2->x[1],&x2->x[2]);
910 flip(&y->x[1],&y->x[2]);
911 flip(&x[1],&x[2]);
912 break;
913 }
914 // one z component is only determined by its radius (without sign)
915 // thus check eight possible sign flips and determine by checking angle with second vector
916 for (i=0;i<8;i++) {
917 // set sign vector accordingly
918 for (j=2;j>=0;j--) {
919 k = (i & pot(2,j)) << j;
920 cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
921 sign[j] = (k == 0) ? 1. : -1.;
922 }
923 cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
924 // apply sign matrix
925 for (j=NDIM;j--;)
926 x[j] *= sign[j];
927 // calculate angle and check
928 ang = x2->Angle (this);
929 cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
930 if (fabs(ang - cos(beta)) < MYEPSILON) {
931 break;
932 }
933 // unapply sign matrix (is its own inverse)
934 for (j=NDIM;j--;)
935 x[j] *= sign[j];
936 }
937 return true;
938};
Note: See TracBrowser for help on using the repository browser.