source: src/molecule_geometry.cpp@ 88104f

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

FIX: CorrelationToSurface() was broken.

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