source: src/molecule_geometry.cpp@ 3839e5

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 3839e5 was 4e10f5, checked in by Tillmann Crueger <crueger@…>, 14 years ago

Merge branch 'stable' into StructureRefactoring

Conflicts:

src/Actions/WorldAction/CenterOnEdgeAction.cpp
src/Actions/WorldAction/ChangeBoxAction.cpp
src/Actions/WorldAction/RepeatBoxAction.cpp
src/Actions/WorldAction/ScaleBoxAction.cpp
src/World.cpp
src/boundary.cpp

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