source: src/molecule_geometry.cpp@ 76096d

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 76096d was 9291d04, checked in by Frederik Heber <heber@…>, 12 years ago

Split saturation with hydrogen of treatment of hydrogen.

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