source: src/molecule_geometry.cpp@ f2cefdb

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 f2cefdb was 112b09, checked in by Tillmann Crueger <crueger@…>, 15 years ago

Added #include "Helpers/MemDebug.hpp" to all .cpp files

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