[e895f7] | 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>
|
---|
[6393ff] | 21 | #include <map>
|
---|
| 22 | #include <set>
|
---|
| 23 | #include <vector>
|
---|
[e895f7] | 24 |
|
---|
| 25 | #include "LinearAlgebra/Vector.hpp"
|
---|
| 26 |
|
---|
[450adf] | 27 | class SphericalPointDistributionTest;
|
---|
| 28 |
|
---|
[e895f7] | 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.) :
|
---|
[0c42f2] | 46 | Bondlength(_Bondlength)
|
---|
[e895f7] | 47 | {}
|
---|
| 48 |
|
---|
| 49 | //!> typedef for the list of points
|
---|
| 50 | typedef std::list<Vector> Polygon_t;
|
---|
[2199c2] | 51 | //!> typedef for the list of points with integral weights
|
---|
| 52 | typedef std::list<std::pair<Vector, int> > WeightedPolygon_t;
|
---|
[6393ff] | 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;
|
---|
[e895f7] | 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 | */
|
---|
[0c42f2] | 64 | template <int N> Polygon_t get() const
|
---|
[e895f7] | 65 | {
|
---|
| 66 | ASSERT(0, "SphericalPointDistribution::get() - not specialized for "+toString(N)+".");
|
---|
[0c42f2] | 67 | return Polygon_t();
|
---|
[e895f7] | 68 | }
|
---|
| 69 |
|
---|
[6393ff] | 70 | template <int N> adjacency_t getConnections()
|
---|
| 71 | {
|
---|
| 72 | ASSERT(0, "SphericalPointDistribution::getConnections() - not specialized for "+toString(N)+".");
|
---|
| 73 | return Polygon_t();
|
---|
| 74 | }
|
---|
| 75 |
|
---|
[3dadd5] | 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 |
|
---|
[fe90ab] | 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.
|
---|
[64cafb2] | 84 | *
|
---|
| 85 | * This is a helper to determine points where to best insert saturation
|
---|
| 86 | * hydrogens.
|
---|
[e895f7] | 87 | *
|
---|
[fe90ab] | 88 | * \param _polygon already filled places to match
|
---|
| 89 | * \param _N desired total number fo points
|
---|
[e895f7] | 90 | */
|
---|
[fe90ab] | 91 | Polygon_t getRemainingPoints(const WeightedPolygon_t &_polygon, const int _N);
|
---|
[64cafb2] | 92 |
|
---|
[e895f7] | 93 | //!> default radius of the spherical distribution
|
---|
| 94 | const double Bondlength;
|
---|
| 95 | //!> precalculated value for root of 3
|
---|
[0c42f2] | 96 | static const double SQRT_3;
|
---|
[653cea] | 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;
|
---|
[b67d89] | 101 |
|
---|
[450adf] | 102 | //!> typedef for a full rotation specification consisting of axis and angle.
|
---|
[b67d89] | 103 | typedef std::pair<Vector, double> Rotation_t;
|
---|
| 104 |
|
---|
[450adf] | 105 | //!> typedef for a list of indices (of points in a polygon)
|
---|
[b67d89] | 106 | typedef std::list<unsigned int> IndexList_t;
|
---|
[450adf] | 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
|
---|
[b67d89] | 110 | typedef std::vector<unsigned int> IndexArray_t;
|
---|
[450adf] | 111 | //!> typedef for a Vector of positions
|
---|
[b67d89] | 112 | typedef std::vector<Vector> VectorArray_t;
|
---|
[450adf] | 113 | //!> typedef for a Vector of positions with weights
|
---|
| 114 | typedef std::vector< std::pair<Vector, int> > WeightedVectorArray_t;
|
---|
[653cea] | 115 | //!> typedef for a vector of degrees (or integral weights)
|
---|
| 116 | typedef std::vector<unsigned int> WeightsArray_t;
|
---|
[b67d89] | 117 |
|
---|
| 118 | //!> amplitude up to which deviations in checks of rotations are tolerated
|
---|
| 119 | static const double warn_amplitude;
|
---|
| 120 |
|
---|
[0983e6] | 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 |
|
---|
[3678eb] | 129 | static Vector calculateCenterOfMinimumDistance(
|
---|
| 130 | const SphericalPointDistribution::VectorArray_t &_positions,
|
---|
| 131 | const SphericalPointDistribution::IndexList_t &_indices);
|
---|
| 132 |
|
---|
[fe90ab] | 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 |
|
---|
[b67d89] | 146 | private:
|
---|
[450adf] | 147 | //!> grant unit tests access to private parts
|
---|
| 148 | friend class SphericalPointDistributionTest;
|
---|
| 149 |
|
---|
[b67d89] | 150 | static std::pair<double, double> calculateErrorOfMatching(
|
---|
[653cea] | 151 | const VectorArray_t &_old,
|
---|
| 152 | const VectorArray_t &_new,
|
---|
| 153 | const IndexTupleList_t &_Matching);
|
---|
[b67d89] | 154 |
|
---|
| 155 | static Polygon_t removeMatchingPoints(
|
---|
[0983e6] | 156 | const PolygonWithIndices &_points);
|
---|
[b67d89] | 157 |
|
---|
| 158 | struct MatchingControlStructure {
|
---|
[fe90ab] | 159 | MatchingControlStructure(
|
---|
| 160 | const adjacency_t &_adjacency,
|
---|
| 161 | const VectorArray_t &_oldpoints,
|
---|
| 162 | const VectorArray_t &_newpoints,
|
---|
| 163 | const WeightsArray_t &_weights
|
---|
| 164 | );
|
---|
[b67d89] | 165 | bool foundflag;
|
---|
| 166 | double bestL2;
|
---|
[fe90ab] | 167 | const adjacency_t &adjacency;
|
---|
| 168 | const VectorArray_t oldpoints;
|
---|
| 169 | const VectorArray_t newpoints;
|
---|
| 170 | const WeightsArray_t weights;
|
---|
[653cea] | 171 | IndexTupleList_t bestmatching;
|
---|
[b67d89] | 172 | };
|
---|
| 173 |
|
---|
| 174 | static void recurseMatchings(
|
---|
| 175 | MatchingControlStructure &_MCS,
|
---|
[653cea] | 176 | IndexTupleList_t &_matching,
|
---|
[b67d89] | 177 | IndexList_t _indices,
|
---|
[653cea] | 178 | WeightsArray_t &_remainingweights,
|
---|
| 179 | WeightsArray_t::iterator _remainiter,
|
---|
| 180 | const unsigned int _matchingsize
|
---|
| 181 | );
|
---|
[b67d89] | 182 |
|
---|
[fe90ab] | 183 | IndexList_t findBestMatching(const WeightedPolygon_t &_polygon);
|
---|
[450adf] | 184 |
|
---|
| 185 | static IndexList_t joinPoints(
|
---|
| 186 | Polygon_t &_newpolygon,
|
---|
| 187 | const VectorArray_t &_newpoints,
|
---|
| 188 | const IndexTupleList_t &_bestmatching
|
---|
[b67d89] | 189 | );
|
---|
| 190 |
|
---|
| 191 | static Rotation_t findPlaneAligningRotation(
|
---|
[0983e6] | 192 | const PolygonWithIndices &_referencepositions,
|
---|
| 193 | const PolygonWithIndices &_currentpositions
|
---|
[b67d89] | 194 | );
|
---|
| 195 |
|
---|
| 196 | static Rotation_t findPointAligningRotation(
|
---|
[0983e6] | 197 | const PolygonWithIndices &remainingold,
|
---|
| 198 | const PolygonWithIndices &remainingnew);
|
---|
[e895f7] | 199 | };
|
---|
| 200 |
|
---|
[f54930] | 201 | // declare specializations
|
---|
| 202 |
|
---|
[0c42f2] | 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;
|
---|
[e895f7] | 217 |
|
---|
[6393ff] | 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 |
|
---|
[e895f7] | 233 | #endif /* SPHERICALPOINTDISTRIBUTION_HPP_ */
|
---|