source: src/molecule_geometry.cpp@ 61ea5b

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 61ea5b was 4e10f5, checked in by Tillmann Crueger <crueger@…>, 15 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
RevLine 
[cee0b57]1/*
2 * molecule_geometry.cpp
3 *
4 * Created on: Oct 5, 2009
5 * Author: heber
6 */
7
[aafd77]8#ifdef HAVE_CONFIG_H
9#include <config.h>
10#endif
11
[112b09]12#include "Helpers/MemDebug.hpp"
13
[f66195]14#include "atom.hpp"
15#include "bond.hpp"
[cee0b57]16#include "config.hpp"
[f66195]17#include "element.hpp"
18#include "helpers.hpp"
19#include "leastsquaremin.hpp"
[36166d]20#include "verbose.hpp"
[e138de]21#include "log.hpp"
[cee0b57]22#include "molecule.hpp"
[b34306]23#include "World.hpp"
[ccf826]24#include "Plane.hpp"
[c94eeb]25#include "Matrix.hpp"
[84c494]26#include "Box.hpp"
[76c0d6]27#include <boost/foreach.hpp>
28
[aafd77]29#include <gsl/gsl_eigen.h>
30#include <gsl/gsl_multimin.h>
31
[cee0b57]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 */
[e138de]39bool molecule::CenterInBox()
[cee0b57]40{
41 bool status = true;
[e138de]42 const Vector *Center = DetermineCenterOfAll();
[eddea2]43 const Vector *CenterBox = DetermineCenterOfBox();
[f429d7]44 Box &domain = World::getInstance().getDomain();
[cee0b57]45
46 // go through all atoms
[273382]47 ActOnAllVectors( &Vector::SubtractVector, *Center);
[eddea2]48 ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
[d0f111]49 BOOST_FOREACH(atom* iter, atoms){
[f429d7]50 *iter->node = domain.WrapPeriodically(*iter->node);
[d0f111]51 }
[cee0b57]52
53 delete(Center);
[52d777]54 delete(CenterBox);
[cee0b57]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 */
[e138de]62bool molecule::BoundInBox()
[cee0b57]63{
64 bool status = true;
[f429d7]65 Box &domain = World::getInstance().getDomain();
[cee0b57]66
67 // go through all atoms
[d0f111]68 BOOST_FOREACH(atom* iter, atoms){
[f429d7]69 *iter->node = domain.WrapPeriodically(*iter->node);
[d0f111]70 }
[cee0b57]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 */
[e138de]79void molecule::CenterEdge(Vector *max)
[cee0b57]80{
81 Vector *min = new Vector;
82
[e138de]83// Log() << Verbose(3) << "Begin of CenterEdge." << endl;
[9879f6]84 molecule::const_iterator iter = begin(); // start at first in list
85 if (iter != end()) { //list not empty?
[cee0b57]86 for (int i=NDIM;i--;) {
[a7b761b]87 max->at(i) = (*iter)->x[i];
88 min->at(i) = (*iter)->x[i];
[cee0b57]89 }
[9879f6]90 for (; iter != end(); ++iter) {// continue with second if present
91 //(*iter)->Output(1,1,out);
[cee0b57]92 for (int i=NDIM;i--;) {
[a7b761b]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);
[cee0b57]95 }
96 }
[e138de]97// Log() << Verbose(4) << "Maximum is ";
[cee0b57]98// max->Output(out);
[e138de]99// Log() << Verbose(0) << ", Minimum is ";
[cee0b57]100// min->Output(out);
[e138de]101// Log() << Verbose(0) << endl;
[cee0b57]102 min->Scale(-1.);
[273382]103 (*max) += (*min);
[cee0b57]104 Translate(min);
105 Center.Zero();
106 }
107 delete(min);
[e138de]108// Log() << Verbose(3) << "End of CenterEdge." << endl;
[cee0b57]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 */
[e138de]115void molecule::CenterOrigin()
[cee0b57]116{
117 int Num = 0;
[9879f6]118 molecule::const_iterator iter = begin(); // start at first in list
[cee0b57]119
120 Center.Zero();
121
[9879f6]122 if (iter != end()) { //list not empty?
123 for (; iter != end(); ++iter) { // continue with second if present
[cee0b57]124 Num++;
[a7b761b]125 Center += (*iter)->x;
[cee0b57]126 }
[bdc91e]127 Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
[cee0b57]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 */
[e138de]136Vector * molecule::DetermineCenterOfAll() const
[cee0b57]137{
[9879f6]138 molecule::const_iterator iter = begin(); // start at first in list
[cee0b57]139 Vector *a = new Vector();
140 double Num = 0;
141
142 a->Zero();
143
[9879f6]144 if (iter != end()) { //list not empty?
145 for (; iter != end(); ++iter) { // continue with second if present
[15b670]146 Num++;
[1024cb]147 (*a) += (*iter)->x;
[cee0b57]148 }
[bdc91e]149 a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
[cee0b57]150 }
151 return a;
152};
153
[eddea2]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);
[84c494]160 const Matrix &M = World::getInstance().getDomain().getM();
[5108e1]161 (*a) *= M;
[eddea2]162 return a;
163};
164
[cee0b57]165/** Returns vector pointing to center of gravity.
166 * \param *out output stream for debugging
167 * \return pointer to center of gravity vector
168 */
[e138de]169Vector * molecule::DetermineCenterOfGravity()
[cee0b57]170{
[9879f6]171 molecule::const_iterator iter = begin(); // start at first in list
[cee0b57]172 Vector *a = new Vector();
173 Vector tmp;
174 double Num = 0;
175
176 a->Zero();
177
[9879f6]178 if (iter != end()) { //list not empty?
179 for (; iter != end(); ++iter) { // continue with second if present
180 Num += (*iter)->type->mass;
[a7b761b]181 tmp = (*iter)->type->mass * (*iter)->x;
[273382]182 (*a) += tmp;
[cee0b57]183 }
[bdc91e]184 a->Scale(1./Num); // divide through total mass
[cee0b57]185 }
[e138de]186// Log() << Verbose(1) << "Resulting center of gravity: ";
[cee0b57]187// a->Output(out);
[e138de]188// Log() << Verbose(0) << endl;
[cee0b57]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 */
[e138de]196void molecule::CenterPeriodic()
[cee0b57]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 */
[e138de]206void molecule::CenterAtVector(Vector *newcenter)
[cee0b57]207{
[273382]208 Center = *newcenter;
[cee0b57]209};
210
211
212/** Scales all atoms by \a *factor.
213 * \param *factor pointer to scaling factor
[1bd79e]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)
[cee0b57]219 */
[776b64]220void molecule::Scale(const double ** const factor)
[cee0b57]221{
[9879f6]222 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[cee0b57]223 for (int j=0;j<MDSteps;j++)
[a7b761b]224 (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
225 (*iter)->x.ScaleAll(*factor);
[cee0b57]226 }
227};
228
229/** Translate all atoms by given vector.
230 * \param trans[] translation vector.
231 */
232void molecule::Translate(const Vector *trans)
233{
[9879f6]234 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[cee0b57]235 for (int j=0;j<MDSteps;j++)
[a7b761b]236 (*iter)->Trajectory.R.at(j) += (*trans);
237 (*iter)->x += (*trans);
[cee0b57]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{
[f429d7]247 Box &domain = World::getInstance().getDomain();
[cee0b57]248
249 // go through all atoms
[eddea2]250 ActOnAllVectors( &Vector::AddVector, *trans);
[d0f111]251 BOOST_FOREACH(atom* iter, atoms){
[f429d7]252 *iter->node = domain.WrapPeriodically(*iter->node);
[d0f111]253 }
[cee0b57]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{
[76c0d6]263 OBSERVE;
[ccf826]264 Plane p(*n,0);
[d0f111]265 BOOST_FOREACH(atom* iter, atoms ){
[76c0d6]266 (*iter->node) = p.mirrorVector(*iter->node);
[ccf826]267 }
[cee0b57]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{
[84c494]275 const Matrix &matrix = World::getInstance().getDomain().getM();
276 const Matrix &inversematrix = World::getInstance().getDomain().getM();
[cee0b57]277 double tmp;
278 bool flag;
279 Vector Testvector, Translationvector;
280
281 do {
282 Center.Zero();
283 flag = true;
[9879f6]284 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[cee0b57]285#ifdef ADDHYDROGEN
[9879f6]286 if ((*iter)->type->Z != 1) {
[cee0b57]287#endif
[5108e1]288 Testvector = inversematrix * (*iter)->x;
[cee0b57]289 Translationvector.Zero();
[9879f6]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
[cee0b57]292 for (int j=0;j<NDIM;j++) {
[a7b761b]293 tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
[cee0b57]294 if ((fabs(tmp)) > BondDistance) {
295 flag = false;
[a7b761b]296 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << (*iter)->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
[cee0b57]297 if (tmp > 0)
[0a4f7f]298 Translationvector[j] -= 1.;
[cee0b57]299 else
[0a4f7f]300 Translationvector[j] += 1.;
[cee0b57]301 }
302 }
303 }
[273382]304 Testvector += Translationvector;
[5108e1]305 Testvector *= matrix;
[273382]306 Center += Testvector;
[0a4f7f]307 Log() << Verbose(1) << "vector is: " << Testvector << endl;
[cee0b57]308#ifdef ADDHYDROGEN
309 // now also change all hydrogens
[9879f6]310 for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
311 if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
[5108e1]312 Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
[273382]313 Testvector += Translationvector;
[5108e1]314 Testvector *= matrix;
[273382]315 Center += Testvector;
[0a4f7f]316 Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
[cee0b57]317 }
318 }
319 }
320#endif
321 }
322 } while (!flag);
[1614174]323
[ea7176]324 Center.Scale(1./static_cast<double>(getAtomCount()));
[cee0b57]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 */
[e138de]332void molecule::PrincipalAxisSystem(bool DoRotate)
[cee0b57]333{
334 double InertiaTensor[NDIM*NDIM];
[e138de]335 Vector *CenterOfGravity = DetermineCenterOfGravity();
[cee0b57]336
[e138de]337 CenterPeriodic();
[cee0b57]338
339 // reset inertia tensor
340 for(int i=0;i<NDIM*NDIM;i++)
341 InertiaTensor[i] = 0.;
342
343 // sum up inertia tensor
[9879f6]344 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[a7b761b]345 Vector x = (*iter)->x;
[cee0b57]346 //x.SubtractVector(CenterOfGravity);
[a7b761b]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]);
[cee0b57]356 }
357 // print InertiaTensor for debugging
[a67d19]358 DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
[cee0b57]359 for(int i=0;i<NDIM;i++) {
360 for(int j=0;j<NDIM;j++)
[a67d19]361 DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
362 DoLog(0) && (Log() << Verbose(0) << endl);
[cee0b57]363 }
[a67d19]364 DoLog(0) && (Log() << Verbose(0) << endl);
[cee0b57]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++) {
[a67d19]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);
[cee0b57]378 }
379
380 // check whether we rotate or not
381 if (DoRotate) {
[a67d19]382 DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
[cee0b57]383 // the eigenvectors specify the transformation matrix
[c94eeb]384 Matrix M = Matrix(evec->data);
[5108e1]385 BOOST_FOREACH(atom* iter, atoms){
386 (*iter->node) *= M;
387 }
[a67d19]388 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
[cee0b57]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
[9879f6]396 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[a7b761b]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]);
[cee0b57]407 }
408 // print InertiaTensor for debugging
[a67d19]409 DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
[cee0b57]410 for(int i=0;i<NDIM;i++) {
411 for(int j=0;j<NDIM;j++)
[a67d19]412 DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
413 DoLog(0) && (Log() << Verbose(0) << endl);
[cee0b57]414 }
[a67d19]415 DoLog(0) && (Log() << Verbose(0) << endl);
[cee0b57]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;
[0a4f7f]432 z_axis[0] = 0.;
433 z_axis[1] = 0.;
434 z_axis[2] = 1.;
[cee0b57]435
436 // rotate on z-x plane
[a67d19]437 DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
[0a4f7f]438 alpha = atan(-n->at(0)/n->at(2));
[a67d19]439 DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
[9879f6]440 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[a7b761b]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];
[cee0b57]444 for (int j=0;j<MDSteps;j++) {
[a7b761b]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];
[cee0b57]448 }
449 }
450 // rotate n vector
[0a4f7f]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);
[8cbb97]454 DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
[cee0b57]455
456 // rotate on z-y plane
[0a4f7f]457 alpha = atan(-n->at(1)/n->at(2));
[a67d19]458 DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
[9879f6]459 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
[a7b761b]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];
[cee0b57]463 for (int j=0;j<MDSteps;j++) {
[a7b761b]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];
[cee0b57]467 }
468 }
469 // rotate n vector (for consistency check)
[0a4f7f]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);
[cee0b57]473
474
[8cbb97]475 DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
[a67d19]476 DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
[cee0b57]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
[0a4f7f]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);
[cee0b57]498 // go through all atoms
[9879f6]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
[a7b761b]501 c = (*iter)->x - a;
[273382]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
[cee0b57]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.