1 /*------------------------------------------------------------------- 2 Copyright 2012 Ravishankar Sundararaman 3 4 This file is part of JDFTx. 5 6 JDFTx is free software: you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 JDFTx is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>. 18 -------------------------------------------------------------------*/ 19 20 #ifndef JDFTX_CORE_WIGNERSEITZ_H 21 #define JDFTX_CORE_WIGNERSEITZ_H 22 23 //! @addtogroup Geometry 24 //! @{ 25 26 #include <core/matrix3.h> 27 #include <float.h> 28 #include <array> 29 #include <list> 30 #include <set> 31 32 //! Wigner-Seitz construction for a 3D lattice (and 2D lattice with orthogonal 3rd direction) 33 class WignerSeitz 34 { 35 public: 36 WignerSeitz(const matrix3<>& R); //!< Construct Wigner-Seitz cell given lattice vectors 37 ~WignerSeitz(); 38 39 //! Find the point within the Wigner-Seitz cell equivalent to x (lattice coordinates) 40 inline vector3<> restrict(const vector3<>& x) const 41 { static const double tol = 1e-8; 42 vector3<> xWS = x; 43 bool changed = true; 44 while(changed) 45 { changed = false; 46 for(const Face* f: faceHalf) 47 { double d = 0.5 * (1. + dot(f->eqn, xWS)); 48 if(d<-tol || d>1.+tol) //not in fundamental zone 49 { xWS -= floor(d) * f->img; 50 changed = true; 51 } 52 } 53 } 54 return xWS; 55 } 56 57 //! Find the point within the Wigner-Seitz cell equivalent to iv (mesh coordinates with sample count S) 58 inline vector3<int> restrict(const vector3<int>& iv, const vector3<int>& S, const vector3<>& invS) const 59 { static const double tol = 1e-8; 60 vector3<int> ivWS = iv; 61 bool changed = true; 62 while(changed) 63 { changed = false; 64 for(const Face* f: faceHalf) 65 { double xDotEqn = 0.; 66 for(int k=0; k<3; k++) 67 xDotEqn += f->eqn[k] * ivWS[k] * invS[k]; 68 double d = 0.5 * (1. + xDotEqn); 69 if(d<-tol || d>1.+tol) //not in fundamental zone 70 { int id = int(floor(d)); 71 for(int k=0; k<3; k++) 72 ivWS[k] -= id * f->img[k] * S[k]; 73 changed = true; 74 } 75 } 76 } 77 return ivWS; 78 } 79 80 //! Find the smallest distance of a point inside the Wigner-Seitz cell from its surface 81 //! Returns 0 if the point is outside the Wigner-Seitz cell 82 //! Ignore direction iDir to obtain 2D behavior, if iDir >= 0 83 inline double boundaryDistance(const vector3<>& x, int iDir=-1) const 84 { double minDistSq = DBL_MAX; 85 for(const Face* f: faceHalf) 86 if(iDir<0 || !f->img[iDir]) 87 { double dDiff = 0.5*(1. - fabs(dot(f->eqn, x))); //fractional distance from planes 88 if(dDiff < 0.) return 0.; //point outside Wigner Seitz cell 89 double distSq = dDiff*dDiff * RTR.metric_length_squared(f->img); 90 if(distSq<minDistSq) minDistSq = distSq; 91 } 92 return sqrt(minDistSq); 93 } 94 95 //! Return true if point is on th eboundray of the Wigner-Seitz cell 96 //! Ignore direction iDir to obtain 2D behavior, if iDir >= 0 97 inline bool onBoundary(const vector3<>& x, int iDir=-1) const 98 { static const double tol = 1e-8; 99 double minDistSq = DBL_MAX; 100 for(const Face* f: faceHalf) 101 if(iDir<0 || !f->img[iDir]) 102 { double dDiff = 0.5*(1. - fabs(dot(f->eqn, x))); //fractional distance from planes 103 if(dDiff < -tol) return false; //point outside Wigner Seitz cell 104 double distSq = dDiff*dDiff * RTR.metric_length_squared(f->img); 105 if(distSq<minDistSq) minDistSq = distSq; 106 } 107 return minDistSq<tol; 108 } 109 110 //! Radius of largest sphere centered at origin contained within the Wigner-Seitz cell 111 //! Ignore direction iDir to obtain 2D behavior, if iDir >= 0 112 double inRadius(int iDir=-1) const; 113 114 //! Radius of smallest sphere centered at origin that contains the Wigner-Seitz cell 115 //! Ignore direction iDir to obtain 2D behavior, if iDir >= 0 116 double circumRadius(int iDir=-1) const; 117 118 //! Get list of neighbouring lattice vectors (along non-truncated directions) that share faces with origin 119 std::vector<vector3<int>> getNeighbours(vector3<bool> isTruncated) const; 120 121 //! Write a wireframe plot to file (for gnuplot) 122 void writeWireframePlot(const char* filename) const; 123 124 //! Write a wireframe plot for Data Explorer (.dx) 125 void writeWireframeDX(const char* filename) const; 126 127 //! Check if the data structure is valid (all links are reciprocated etc.) 128 void checkGraph() const; 129 130 //! Output vertex, edge and face connectivity info: 131 void writeGraph(FILE* fp=stdout) const; 132 133 static const double geomRelTol; //!< relative tolerance for orthogonality and volume checks 134 135 //! Check whether lattice vectors are orthogonal (within relative tolerance geomRelTol) 136 static bool isOrthogonal(const vector3<>&, const vector3<>&); 137 138 private: 139 struct Vertex; //!< Point 140 struct Edge; //!< Line segment 141 struct Face; //!<Polygonal facet 142 143 matrix3<> R, invR, RTR; //!< matrix of lattice vectors, and its combinations 144 double minDistSq; //!< threshold on distance squared for welding vertices 145 std::list<Vertex*> vertex; //!< set of all vertices 146 std::set<Edge*> edge; //!< set of all edges 147 std::set<Face*> face; //!< set of all faces 148 std::vector<Face*> faceHalf; //!< array of half the faces, one from each inversion symmetry pair 149 150 //! Point 151 struct Vertex 152 { vector3<> pos; //!< position in lattice coordinates 153 std::list<Edge*> edge; //!< edges bounded by this vertex 154 }; 155 156 //! Line segment 157 struct Edge 158 { std::array<Vertex*,2> vertex; //!< vertices bounding this edge 159 std::array<Face*,2> face; //!< faces bounded by this edge 160 }; 161 162 //! Polygonal facet 163 struct Face 164 { vector3<int> img; //!< image of origin under the plane containing this face (lattice coordinates) 165 vector3<> eqn; //!< equation of plane given by eqn.x==1 (x in lattice coordinates) 166 std::list<Edge*> edge; //!< edges bounding this face 167 }; 168 169 //! Slice the current polyhedron by the perpendicular bisector of 0->a (a in lattice coordinates) 170 void addPlane(const vector3<int>& a); 171 172 //! Add a edge from vStart towards vEnd in face f 173 //! Note: edges are added to the end of the face list by default (therefore must call in order) 174 //! However if or if there is only one missing edge (lastEdge=true), insertion will be at correct location 175 void addEdge(Face* f, Vertex* vStart, Vertex* vEnd, bool lastEdge=false); 176 }; 177 178 //! @} 179 #endif // JDFTX_CORE_WIGNERSEITZ_H 180