source: src/molecule_geometry.cpp@ 8453b3

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 8453b3 was f221e3, checked in by Frederik Heber <heber@…>, 13 years ago

FIX: molecule::RotateToPrincipalAxisSystem() now always normalizes Axis vector.

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