1 | /*
|
---|
2 | * SphericalPointDistribution.hpp
|
---|
3 | *
|
---|
4 | * Created on: May 29, 2014
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 |
|
---|
9 | #ifndef SPHERICALPOINTDISTRIBUTION_HPP_
|
---|
10 | #define SPHERICALPOINTDISTRIBUTION_HPP_
|
---|
11 |
|
---|
12 | // include config.h
|
---|
13 | #ifdef HAVE_CONFIG_H
|
---|
14 | #include <config.h>
|
---|
15 | #endif
|
---|
16 |
|
---|
17 | #include "CodePatterns/Assert.hpp"
|
---|
18 |
|
---|
19 | #include <cmath>
|
---|
20 | #include <list>
|
---|
21 | #include <map>
|
---|
22 | #include <set>
|
---|
23 | #include <vector>
|
---|
24 |
|
---|
25 | #include "LinearAlgebra/Vector.hpp"
|
---|
26 |
|
---|
27 | class SphericalPointDistributionTest;
|
---|
28 |
|
---|
29 | /** contains getters for the VSEPR model for specific number of electrons.
|
---|
30 | *
|
---|
31 | * This struct contains specialized functions returning a list of Vectors
|
---|
32 | * (points in space) to match the VSEPR model for the given number of electrons.
|
---|
33 | *
|
---|
34 | * This is implemented via template specialization of the function get().
|
---|
35 | *
|
---|
36 | * These specializations are taken from the python script \b CreateVspeShapes.py
|
---|
37 | * by Christian Neuen, 07th May 2009.
|
---|
38 | */
|
---|
39 | struct SphericalPointDistribution
|
---|
40 | {
|
---|
41 | /** Cstor for SphericalPointDistribution, allows setting radius of sphere
|
---|
42 | *
|
---|
43 | * \param _BondLength desired radius of sphere
|
---|
44 | */
|
---|
45 | SphericalPointDistribution(const double _Bondlength = 1.) :
|
---|
46 | Bondlength(_Bondlength)
|
---|
47 | {}
|
---|
48 |
|
---|
49 | //!> typedef for the list of points
|
---|
50 | typedef std::list<Vector> Polygon_t;
|
---|
51 | //!> typedef for the list of points with integral weights
|
---|
52 | typedef std::list<std::pair<Vector, int> > WeightedPolygon_t;
|
---|
53 | //!> typedef for a sorted list of indices
|
---|
54 | typedef std::set<unsigned int> IndexSet_t;
|
---|
55 | //!> typedef for the adjacency list of a polygon
|
---|
56 | typedef std::map<unsigned int, IndexSet_t > adjacency_t;
|
---|
57 |
|
---|
58 | /** General getter function for the distribution of points on the surface.
|
---|
59 | *
|
---|
60 | * \warn this function needs to be specialized!
|
---|
61 | *
|
---|
62 | * \return Polygon_t with points on the surface centered at (0,0,0)
|
---|
63 | */
|
---|
64 | template <int N> Polygon_t get() const
|
---|
65 | {
|
---|
66 | ASSERT(0, "SphericalPointDistribution::get() - not specialized for "+toString(N)+".");
|
---|
67 | return Polygon_t();
|
---|
68 | }
|
---|
69 |
|
---|
70 | template <int N> adjacency_t getConnections()
|
---|
71 | {
|
---|
72 | ASSERT(0, "SphericalPointDistribution::getConnections() - not specialized for "+toString(N)+".");
|
---|
73 | return Polygon_t();
|
---|
74 | }
|
---|
75 |
|
---|
76 | /** Initializes the polygon with the given \a _NumberOfPoints.
|
---|
77 | *
|
---|
78 | * \param _NumberOfPoints number of points
|
---|
79 | */
|
---|
80 | Polygon_t getSimplePolygon(const int _NumberOfPoints) const;
|
---|
81 |
|
---|
82 | /** Returns vacant spots to fill to get a complete spherical point distribution from
|
---|
83 | * given points \a _polygon, containing then \a _N in total.
|
---|
84 | *
|
---|
85 | * This is a helper to determine points where to best insert saturation
|
---|
86 | * hydrogens.
|
---|
87 | *
|
---|
88 | * \param _polygon already filled places to match
|
---|
89 | * \param _N desired total number fo points
|
---|
90 | */
|
---|
91 | Polygon_t getRemainingPoints(const WeightedPolygon_t &_polygon, const int _N);
|
---|
92 |
|
---|
93 | //!> default radius of the spherical distribution
|
---|
94 | const double Bondlength;
|
---|
95 | //!> precalculated value for root of 3
|
---|
96 | static const double SQRT_3;
|
---|
97 | //!> threshold for L1 error below which matching is immediately acceptable
|
---|
98 | static const double L1THRESHOLD;
|
---|
99 | //!> threshold for L2 error below which matching is acceptable
|
---|
100 | static const double L2THRESHOLD;
|
---|
101 |
|
---|
102 | //!> typedef for a full rotation specification consisting of axis and angle.
|
---|
103 | typedef std::pair<Vector, double> Rotation_t;
|
---|
104 |
|
---|
105 | //!> typedef for a list of indices (of points in a polygon)
|
---|
106 | typedef std::list<unsigned int> IndexList_t;
|
---|
107 | //!> typedef enumerating possibly multiple points accumulated as one point
|
---|
108 | typedef std::list< IndexList_t > IndexTupleList_t;
|
---|
109 | //!> typedef for a vector of indices
|
---|
110 | typedef std::vector<unsigned int> IndexArray_t;
|
---|
111 | //!> typedef for a Vector of positions
|
---|
112 | typedef std::vector<Vector> VectorArray_t;
|
---|
113 | //!> typedef for a Vector of positions with weights
|
---|
114 | typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
|
---|
115 | //!> typedef for a vector of degrees (or integral weights)
|
---|
116 | typedef std::vector<unsigned int> WeightsArray_t;
|
---|
117 |
|
---|
118 | //!> amplitude up to which deviations in checks of rotations are tolerated
|
---|
119 | static const double warn_amplitude;
|
---|
120 |
|
---|
121 | struct PolygonWithIndices
|
---|
122 | {
|
---|
123 | //!> array with points
|
---|
124 | VectorArray_t polygon;
|
---|
125 | //!> list with indices for the above points, defining subset
|
---|
126 | IndexList_t indices;
|
---|
127 | };
|
---|
128 |
|
---|
129 | static Vector calculateCenterOfMinimumDistance(
|
---|
130 | const SphericalPointDistribution::VectorArray_t &_positions,
|
---|
131 | const SphericalPointDistribution::IndexList_t &_indices);
|
---|
132 |
|
---|
133 | private:
|
---|
134 | //!> points for the ideal distribution
|
---|
135 | Polygon_t points;
|
---|
136 | //!> connection information between these ideal points
|
---|
137 | adjacency_t adjacency;
|
---|
138 |
|
---|
139 | /** Initialize inner status (points and adjacency) to desired number of
|
---|
140 | * points.
|
---|
141 | *
|
---|
142 | * \param _N number of points
|
---|
143 | */
|
---|
144 | void initSelf(const int _N);
|
---|
145 |
|
---|
146 | private:
|
---|
147 | //!> grant unit tests access to private parts
|
---|
148 | friend class SphericalPointDistributionTest;
|
---|
149 |
|
---|
150 | static std::pair<double, double> calculateErrorOfMatching(
|
---|
151 | const VectorArray_t &_old,
|
---|
152 | const VectorArray_t &_new,
|
---|
153 | const IndexTupleList_t &_Matching);
|
---|
154 |
|
---|
155 | static Polygon_t removeMatchingPoints(
|
---|
156 | const PolygonWithIndices &_points);
|
---|
157 |
|
---|
158 | struct MatchingControlStructure {
|
---|
159 | MatchingControlStructure(
|
---|
160 | const adjacency_t &_adjacency,
|
---|
161 | const VectorArray_t &_oldpoints,
|
---|
162 | const VectorArray_t &_newpoints,
|
---|
163 | const WeightsArray_t &_weights
|
---|
164 | );
|
---|
165 | bool foundflag;
|
---|
166 | double bestL2;
|
---|
167 | const adjacency_t &adjacency;
|
---|
168 | const VectorArray_t oldpoints;
|
---|
169 | const VectorArray_t newpoints;
|
---|
170 | const WeightsArray_t weights;
|
---|
171 | IndexTupleList_t bestmatching;
|
---|
172 | };
|
---|
173 |
|
---|
174 | static void recurseMatchings(
|
---|
175 | MatchingControlStructure &_MCS,
|
---|
176 | IndexTupleList_t &_matching,
|
---|
177 | IndexList_t _indices,
|
---|
178 | WeightsArray_t &_remainingweights,
|
---|
179 | WeightsArray_t::iterator _remainiter,
|
---|
180 | const unsigned int _matchingsize
|
---|
181 | );
|
---|
182 |
|
---|
183 | IndexList_t findBestMatching(const WeightedPolygon_t &_polygon);
|
---|
184 |
|
---|
185 | static IndexList_t joinPoints(
|
---|
186 | Polygon_t &_newpolygon,
|
---|
187 | const VectorArray_t &_newpoints,
|
---|
188 | const IndexTupleList_t &_bestmatching
|
---|
189 | );
|
---|
190 |
|
---|
191 | static Rotation_t findPlaneAligningRotation(
|
---|
192 | const PolygonWithIndices &_referencepositions,
|
---|
193 | const PolygonWithIndices &_currentpositions
|
---|
194 | );
|
---|
195 |
|
---|
196 | static Rotation_t findPointAligningRotation(
|
---|
197 | const PolygonWithIndices &remainingold,
|
---|
198 | const PolygonWithIndices &remainingnew);
|
---|
199 | };
|
---|
200 |
|
---|
201 | // declare specializations
|
---|
202 |
|
---|
203 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>() const;
|
---|
204 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>() const;
|
---|
205 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>() const;
|
---|
206 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>() const;
|
---|
207 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>() const;
|
---|
208 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>() const;
|
---|
209 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>() const;
|
---|
210 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>() const;
|
---|
211 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>() const;
|
---|
212 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>() const;
|
---|
213 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>() const;
|
---|
214 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>() const;
|
---|
215 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>() const;
|
---|
216 | template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>() const;
|
---|
217 |
|
---|
218 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<0>();
|
---|
219 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<1>();
|
---|
220 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<2>();
|
---|
221 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<3>();
|
---|
222 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<4>();
|
---|
223 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<5>();
|
---|
224 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<6>();
|
---|
225 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<7>();
|
---|
226 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<8>();
|
---|
227 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<9>();
|
---|
228 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<10>();
|
---|
229 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<11>();
|
---|
230 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<12>();
|
---|
231 | template <> SphericalPointDistribution::adjacency_t SphericalPointDistribution::getConnections<14>();
|
---|
232 |
|
---|
233 | #endif /* SPHERICALPOINTDISTRIBUTION_HPP_ */
|
---|