source: src/molecule_geometry.cpp@ eddea2

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 eddea2 was eddea2, checked in by Frederik Heber <heber@…>, 15 years ago

Added case '-T' (periodic translation) to testsuite.

  • BUGFIX: Vector::WrapPeriodically() added half of the box on top which is not as intended.
    • this broke testcase '-b', which relied on the half of the box adding.
    • new function molecule::DetermineCenterOfBox()
    • is used instead to additionally add this center in molecule::CenterInBox()
  • BUGFIX: molecule::TranslatePeriodically() subtracted vector instead of adding.
  • Property mode set to 100644
File size: 18.0 KB
Line 
1/*
2 * molecule_geometry.cpp
3 *
4 * Created on: Oct 5, 2009
5 * Author: heber
6 */
7
8#include "atom.hpp"
9#include "bond.hpp"
10#include "config.hpp"
11#include "element.hpp"
12#include "helpers.hpp"
13#include "leastsquaremin.hpp"
14#include "log.hpp"
15#include "memoryallocator.hpp"
16#include "molecule.hpp"
17#include "World.hpp"
18
19/************************************* Functions for class molecule *********************************/
20
21
22/** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths.
23 * \param *out output stream for debugging
24 */
25bool molecule::CenterInBox()
26{
27 bool status = true;
28 const Vector *Center = DetermineCenterOfAll();
29 const Vector *CenterBox = DetermineCenterOfBox();
30 double * const cell_size = World::getInstance().getDomain();
31 double *M = ReturnFullMatrixforSymmetric(cell_size);
32 double *Minv = InverseMatrix(M);
33
34 // go through all atoms
35 ActOnAllVectors( &Vector::SubtractVector, *Center);
36 ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
37 ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
38
39 Free(&M);
40 Free(&Minv);
41 delete(Center);
42 return status;
43};
44
45
46/** Bounds the molecule in the box whose lengths are defined by vector \a *BoxLengths.
47 * \param *out output stream for debugging
48 */
49bool molecule::BoundInBox()
50{
51 bool status = true;
52 double * const cell_size = World::getInstance().getDomain();
53 double *M = ReturnFullMatrixforSymmetric(cell_size);
54 double *Minv = InverseMatrix(M);
55
56 // go through all atoms
57 ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
58
59 Free(&M);
60 Free(&Minv);
61 return status;
62};
63
64/** Centers the edge of the atoms at (0,0,0).
65 * \param *out output stream for debugging
66 * \param *max coordinates of other edge, specifying box dimensions.
67 */
68void molecule::CenterEdge(Vector *max)
69{
70 Vector *min = new Vector;
71
72// Log() << Verbose(3) << "Begin of CenterEdge." << endl;
73 atom *ptr = start->next; // start at first in list
74 if (ptr != end) { //list not empty?
75 for (int i=NDIM;i--;) {
76 max->at(i) = ptr->x[i];
77 min->at(i) = ptr->x[i];
78 }
79 while (ptr->next != end) { // continue with second if present
80 ptr = ptr->next;
81 //ptr->Output(1,1,out);
82 for (int i=NDIM;i--;) {
83 max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
84 min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(i);
85 }
86 }
87// Log() << Verbose(4) << "Maximum is ";
88// max->Output(out);
89// Log() << Verbose(0) << ", Minimum is ";
90// min->Output(out);
91// Log() << Verbose(0) << endl;
92 min->Scale(-1.);
93 (*max) += (*min);
94 Translate(min);
95 Center.Zero();
96 }
97 delete(min);
98// Log() << Verbose(3) << "End of CenterEdge." << endl;
99};
100
101/** Centers the center of the atoms at (0,0,0).
102 * \param *out output stream for debugging
103 * \param *center return vector for translation vector
104 */
105void molecule::CenterOrigin()
106{
107 int Num = 0;
108 atom *ptr = start; // start at first in list
109
110 Center.Zero();
111
112 if (ptr->next != end) { //list not empty?
113 while (ptr->next != end) {
114 ptr = ptr->next;
115 Num++;
116 Center += ptr->x;
117 }
118 Center.Scale(-1./Num); // divide through total number (and sign for direction)
119 Translate(&Center);
120 Center.Zero();
121 }
122};
123
124/** Returns vector pointing to center of all atoms.
125 * \return pointer to center of all vector
126 */
127Vector * molecule::DetermineCenterOfAll() const
128{
129 atom *ptr = start; // start at first in list
130 Vector *a = new Vector();
131 double Num = 0;
132
133 a->Zero();
134
135 if (ptr->next != end) { //list not empty?
136 while (ptr->next != end) {
137 ptr = ptr->next;
138 Num += 1.;
139 (*a) += ptr->x;
140 }
141 a->Scale(1./Num); // divide through total mass (and sign for direction)
142 }
143 return a;
144};
145
146/** Returns vector pointing to center of the domain.
147 * \return pointer to center of the domain
148 */
149Vector * molecule::DetermineCenterOfBox() const
150{
151 Vector *a = new Vector(0.5,0.5,0.5);
152
153 const double *cell_size = World::getInstance().getDomain();
154 double *M = ReturnFullMatrixforSymmetric(cell_size);
155 a->MatrixMultiplication(M);
156 Free(&M);
157
158 return a;
159};
160
161/** Returns vector pointing to center of gravity.
162 * \param *out output stream for debugging
163 * \return pointer to center of gravity vector
164 */
165Vector * molecule::DetermineCenterOfGravity()
166{
167 atom *ptr = start->next; // start at first in list
168 Vector *a = new Vector();
169 Vector tmp;
170 double Num = 0;
171
172 a->Zero();
173
174 if (ptr != end) { //list not empty?
175 while (ptr->next != end) { // continue with second if present
176 ptr = ptr->next;
177 Num += ptr->type->mass;
178 tmp = ptr->type->mass * ptr->x;
179 (*a) += tmp;
180 }
181 a->Scale(-1./Num); // divide through total mass (and sign for direction)
182 }
183// Log() << Verbose(1) << "Resulting center of gravity: ";
184// a->Output(out);
185// Log() << Verbose(0) << endl;
186 return a;
187};
188
189/** Centers the center of gravity of the atoms at (0,0,0).
190 * \param *out output stream for debugging
191 * \param *center return vector for translation vector
192 */
193void molecule::CenterPeriodic()
194{
195 DeterminePeriodicCenter(Center);
196};
197
198
199/** Centers the center of gravity of the atoms at (0,0,0).
200 * \param *out output stream for debugging
201 * \param *center return vector for translation vector
202 */
203void molecule::CenterAtVector(Vector *newcenter)
204{
205 Center = *newcenter;
206};
207
208
209/** Scales all atoms by \a *factor.
210 * \param *factor pointer to scaling factor
211 *
212 * TODO: Is this realy what is meant, i.e.
213 * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
214 * or rather
215 * x=(**factor) * x (as suggested by comment)
216 */
217void molecule::Scale(const double ** const factor)
218{
219 atom *ptr = start;
220
221 while (ptr->next != end) {
222 ptr = ptr->next;
223 for (int j=0;j<MDSteps;j++)
224 ptr->Trajectory.R.at(j).ScaleAll(*factor);
225 ptr->x.ScaleAll(*factor);
226 }
227};
228
229/** Translate all atoms by given vector.
230 * \param trans[] translation vector.
231 */
232void molecule::Translate(const Vector *trans)
233{
234 atom *ptr = start;
235
236 while (ptr->next != end) {
237 ptr = ptr->next;
238 for (int j=0;j<MDSteps;j++)
239 ptr->Trajectory.R.at(j) += (*trans);
240 ptr->x += (*trans);
241 }
242};
243
244/** Translate the molecule periodically in the box.
245 * \param trans[] translation vector.
246 * TODO treatment of trajetories missing
247 */
248void molecule::TranslatePeriodically(const Vector *trans)
249{
250 double * const cell_size = World::getInstance().getDomain();
251 double *M = ReturnFullMatrixforSymmetric(cell_size);
252 double *Minv = InverseMatrix(M);
253
254 // go through all atoms
255 ActOnAllVectors( &Vector::AddVector, *trans);
256 ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
257
258 Free(&M);
259 Free(&Minv);
260};
261
262
263/** Mirrors all atoms against a given plane.
264 * \param n[] normal vector of mirror plane.
265 */
266void molecule::Mirror(const Vector *n)
267{
268 ActOnAllVectors( &Vector::Mirror, *n );
269};
270
271/** Determines center of molecule (yet not considering atom masses).
272 * \param center reference to return vector
273 */
274void molecule::DeterminePeriodicCenter(Vector &center)
275{
276 atom *Walker = start;
277 double * const cell_size = World::getInstance().getDomain();
278 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
279 double *inversematrix = InverseMatrix(cell_size);
280 double tmp;
281 bool flag;
282 Vector Testvector, Translationvector;
283
284 do {
285 Center.Zero();
286 flag = true;
287 while (Walker->next != end) {
288 Walker = Walker->next;
289#ifdef ADDHYDROGEN
290 if (Walker->type->Z != 1) {
291#endif
292 Testvector = Walker->x;
293 Testvector.MatrixMultiplication(inversematrix);
294 Translationvector.Zero();
295 for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
296 if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
297 for (int j=0;j<NDIM;j++) {
298 tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
299 if ((fabs(tmp)) > BondDistance) {
300 flag = false;
301 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
302 if (tmp > 0)
303 Translationvector[j] -= 1.;
304 else
305 Translationvector[j] += 1.;
306 }
307 }
308 }
309 Testvector += Translationvector;
310 Testvector.MatrixMultiplication(matrix);
311 Center += Testvector;
312 Log() << Verbose(1) << "vector is: " << Testvector << endl;
313#ifdef ADDHYDROGEN
314 // now also change all hydrogens
315 for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
316 if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
317 Testvector = (*Runner)->GetOtherAtom(Walker)->x;
318 Testvector.MatrixMultiplication(inversematrix);
319 Testvector += Translationvector;
320 Testvector.MatrixMultiplication(matrix);
321 Center += Testvector;
322 Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
323 }
324 }
325 }
326#endif
327 }
328 } while (!flag);
329 Free(&matrix);
330 Free(&inversematrix);
331
332 Center.Scale(1./(double)AtomCount);
333};
334
335/** Transforms/Rotates the given molecule into its principal axis system.
336 * \param *out output stream for debugging
337 * \param DoRotate whether to rotate (true) or only to determine the PAS.
338 * TODO treatment of trajetories missing
339 */
340void molecule::PrincipalAxisSystem(bool DoRotate)
341{
342 atom *ptr = start; // start at first in list
343 double InertiaTensor[NDIM*NDIM];
344 Vector *CenterOfGravity = DetermineCenterOfGravity();
345
346 CenterPeriodic();
347
348 // reset inertia tensor
349 for(int i=0;i<NDIM*NDIM;i++)
350 InertiaTensor[i] = 0.;
351
352 // sum up inertia tensor
353 while (ptr->next != end) {
354 ptr = ptr->next;
355 Vector x = ptr->x;
356 //x.SubtractVector(CenterOfGravity);
357 InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
358 InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
359 InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
360 InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
361 InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
362 InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
363 InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
364 InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
365 InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
366 }
367 // print InertiaTensor for debugging
368 DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
369 for(int i=0;i<NDIM;i++) {
370 for(int j=0;j<NDIM;j++)
371 DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
372 DoLog(0) && (Log() << Verbose(0) << endl);
373 }
374 DoLog(0) && (Log() << Verbose(0) << endl);
375
376 // diagonalize to determine principal axis system
377 gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
378 gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
379 gsl_vector *eval = gsl_vector_alloc(NDIM);
380 gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
381 gsl_eigen_symmv(&m.matrix, eval, evec, T);
382 gsl_eigen_symmv_free(T);
383 gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
384
385 for(int i=0;i<NDIM;i++) {
386 DoLog(1) && (Log() << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i));
387 DoLog(0) && (Log() << Verbose(0) << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl);
388 }
389
390 // check whether we rotate or not
391 if (DoRotate) {
392 DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
393 // the eigenvectors specify the transformation matrix
394 ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
395 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
396
397 // summing anew for debugging (resulting matrix has to be diagonal!)
398 // reset inertia tensor
399 for(int i=0;i<NDIM*NDIM;i++)
400 InertiaTensor[i] = 0.;
401
402 // sum up inertia tensor
403 ptr = start;
404 while (ptr->next != end) {
405 ptr = ptr->next;
406 Vector x = ptr->x;
407 //x.SubtractVector(CenterOfGravity);
408 InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
409 InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
410 InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
411 InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
412 InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
413 InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
414 InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
415 InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
416 InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
417 }
418 // print InertiaTensor for debugging
419 DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
420 for(int i=0;i<NDIM;i++) {
421 for(int j=0;j<NDIM;j++)
422 DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
423 DoLog(0) && (Log() << Verbose(0) << endl);
424 }
425 DoLog(0) && (Log() << Verbose(0) << endl);
426 }
427
428 // free everything
429 delete(CenterOfGravity);
430 gsl_vector_free(eval);
431 gsl_matrix_free(evec);
432};
433
434
435/** Align all atoms in such a manner that given vector \a *n is along z axis.
436 * \param n[] alignment vector.
437 */
438void molecule::Align(Vector *n)
439{
440 atom *ptr = start;
441 double alpha, tmp;
442 Vector z_axis;
443 z_axis[0] = 0.;
444 z_axis[1] = 0.;
445 z_axis[2] = 1.;
446
447 // rotate on z-x plane
448 DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
449 alpha = atan(-n->at(0)/n->at(2));
450 DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
451 while (ptr->next != end) {
452 ptr = ptr->next;
453 tmp = ptr->x[0];
454 ptr->x[0] = cos(alpha) * tmp + sin(alpha) * ptr->x[2];
455 ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
456 for (int j=0;j<MDSteps;j++) {
457 tmp = ptr->Trajectory.R.at(j)[0];
458 ptr->Trajectory.R.at(j)[0] = cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
459 ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
460 }
461 }
462 // rotate n vector
463 tmp = n->at(0);
464 n->at(0) = cos(alpha) * tmp + sin(alpha) * n->at(2);
465 n->at(2) = -sin(alpha) * tmp + cos(alpha) * n->at(2);
466 DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
467
468 // rotate on z-y plane
469 ptr = start;
470 alpha = atan(-n->at(1)/n->at(2));
471 DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
472 while (ptr->next != end) {
473 ptr = ptr->next;
474 tmp = ptr->x[1];
475 ptr->x[1] = cos(alpha) * tmp + sin(alpha) * ptr->x[2];
476 ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
477 for (int j=0;j<MDSteps;j++) {
478 tmp = ptr->Trajectory.R.at(j)[1];
479 ptr->Trajectory.R.at(j)[1] = cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
480 ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
481 }
482 }
483 // rotate n vector (for consistency check)
484 tmp = n->at(1);
485 n->at(1) = cos(alpha) * tmp + sin(alpha) * n->at(2);
486 n->at(2) = -sin(alpha) * tmp + cos(alpha) * n->at(2);
487
488
489 DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
490 DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
491};
492
493
494/** Calculates sum over least square distance to line hidden in \a *x.
495 * \param *x offset and direction vector
496 * \param *params pointer to lsq_params structure
497 * \return \f$ sum_i^N | y_i - (a + t_i b)|^2\f$
498 */
499double LeastSquareDistance (const gsl_vector * x, void * params)
500{
501 double res = 0, t;
502 Vector a,b,c,d;
503 struct lsq_params *par = (struct lsq_params *)params;
504 atom *ptr = par->mol->start;
505
506 // initialize vectors
507 a[0] = gsl_vector_get(x,0);
508 a[1] = gsl_vector_get(x,1);
509 a[2] = gsl_vector_get(x,2);
510 b[0] = gsl_vector_get(x,3);
511 b[1] = gsl_vector_get(x,4);
512 b[2] = gsl_vector_get(x,5);
513 // go through all atoms
514 while (ptr != par->mol->end) {
515 ptr = ptr->next;
516 if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
517 c = ptr->x - a;
518 t = c.ScalarProduct(b); // get direction parameter
519 d = t*b; // and create vector
520 c -= d; // ... yielding distance vector
521 res += d.ScalarProduct(d); // add squared distance
522 }
523 }
524 return res;
525};
526
527/** By minimizing the least square distance gains alignment vector.
528 * \bug this is not yet working properly it seems
529 */
530void molecule::GetAlignvector(struct lsq_params * par) const
531{
532 int np = 6;
533
534 const gsl_multimin_fminimizer_type *T =
535 gsl_multimin_fminimizer_nmsimplex;
536 gsl_multimin_fminimizer *s = NULL;
537 gsl_vector *ss;
538 gsl_multimin_function minex_func;
539
540 size_t iter = 0, i;
541 int status;
542 double size;
543
544 /* Initial vertex size vector */
545 ss = gsl_vector_alloc (np);
546
547 /* Set all step sizes to 1 */
548 gsl_vector_set_all (ss, 1.0);
549
550 /* Starting point */
551 par->x = gsl_vector_alloc (np);
552 par->mol = this;
553
554 gsl_vector_set (par->x, 0, 0.0); // offset
555 gsl_vector_set (par->x, 1, 0.0);
556 gsl_vector_set (par->x, 2, 0.0);
557 gsl_vector_set (par->x, 3, 0.0); // direction
558 gsl_vector_set (par->x, 4, 0.0);
559 gsl_vector_set (par->x, 5, 1.0);
560
561 /* Initialize method and iterate */
562 minex_func.f = &LeastSquareDistance;
563 minex_func.n = np;
564 minex_func.params = (void *)par;
565
566 s = gsl_multimin_fminimizer_alloc (T, np);
567 gsl_multimin_fminimizer_set (s, &minex_func, par->x, ss);
568
569 do
570 {
571 iter++;
572 status = gsl_multimin_fminimizer_iterate(s);
573
574 if (status)
575 break;
576
577 size = gsl_multimin_fminimizer_size (s);
578 status = gsl_multimin_test_size (size, 1e-2);
579
580 if (status == GSL_SUCCESS)
581 {
582 printf ("converged to minimum at\n");
583 }
584
585 printf ("%5d ", (int)iter);
586 for (i = 0; i < (size_t)np; i++)
587 {
588 printf ("%10.3e ", gsl_vector_get (s->x, i));
589 }
590 printf ("f() = %7.3f size = %.3f\n", s->fval, size);
591 }
592 while (status == GSL_CONTINUE && iter < 100);
593
594 for (i=0;i<(size_t)np;i++)
595 gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
596 //gsl_vector_free(par->x);
597 gsl_vector_free(ss);
598 gsl_multimin_fminimizer_free (s);
599};
Note: See TracBrowser for help on using the repository browser.