WignerSeitz.h

Go to the documentation of this file.

/*-------------------------------------------------------------------
Copyright 2012 Ravishankar Sundararaman

This file is part of JDFTx.

JDFTx is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

JDFTx is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with JDFTx.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------*/

#ifndef JDFTX_CORE_WIGNERSEITZ_H
#define JDFTX_CORE_WIGNERSEITZ_H


#include <core/matrix3.h>
#include <float.h>
#include <array>
#include <list>
#include <set>

class WignerSeitz
{
public:
        WignerSeitz(const matrix3<>& R); 
        ~WignerSeitz();

        inline vector3<> restrict(const vector3<>& x) const
        {       static const double tol = 1e-8;
                vector3<> xWS = x;
                bool changed = true;
                while(changed)
                {       changed = false;
                        for(const Face* f: faceHalf)
                        {       double d = 0.5 * (1. + dot(f->eqn, xWS));
                                if(d<-tol || d>1.+tol) //not in fundamental zone
                                {       xWS -= floor(d) * f->img;
                                        changed = true;
                                }
                        }
                }
                return xWS;
        }
        
        inline vector3<int> restrict(const vector3<int>& iv, const vector3<int>& S, const vector3<>& invS) const
        {       static const double tol = 1e-8;
                vector3<int> ivWS = iv;
                bool changed = true;
                while(changed)
                {       changed = false;
                        for(const Face* f: faceHalf)
                        {       double xDotEqn = 0.;
                                for(int k=0; k<3; k++)
                                        xDotEqn += f->eqn[k] * ivWS[k] * invS[k];
                                double d = 0.5 * (1. + xDotEqn);
                                if(d<-tol || d>1.+tol) //not in fundamental zone
                                {       int id = int(floor(d));
                                        for(int k=0; k<3; k++)
                                                ivWS[k] -= id * f->img[k] * S[k];
                                        changed = true;
                                }
                        }
                }
                return ivWS;
        }
        
        inline double boundaryDistance(const vector3<>& x, int iDir=-1) const
        {       double minDistSq = DBL_MAX;
                for(const Face* f: faceHalf)
                        if(iDir<0 || !f->img[iDir])
                        {       double dDiff = 0.5*(1. - fabs(dot(f->eqn, x))); //fractional distance from planes
                                if(dDiff < 0.) return 0.; //point outside Wigner Seitz cell
                                double distSq = dDiff*dDiff * RTR.metric_length_squared(f->img);
                                if(distSq<minDistSq) minDistSq = distSq;
                        }
                return sqrt(minDistSq);
        }
        
        inline bool onBoundary(const vector3<>& x, int iDir=-1) const
        {       static const double tol = 1e-8;
                double minDistSq = DBL_MAX;
                for(const Face* f: faceHalf)
                        if(iDir<0 || !f->img[iDir])
                        {       double dDiff = 0.5*(1. - fabs(dot(f->eqn, x))); //fractional distance from planes
                                if(dDiff < -tol) return false; //point outside Wigner Seitz cell
                                double distSq = dDiff*dDiff * RTR.metric_length_squared(f->img);
                                if(distSq<minDistSq) minDistSq = distSq;
                        }
                return minDistSq<tol;
        }

        double inRadius(int iDir=-1) const;
        
        double circumRadius(int iDir=-1) const;
        
        void writeWireframePlot(const char* filename) const;
        
        void writeWireframeDX(const char* filename) const;

        void checkGraph() const;
        
        void writeGraph(FILE* fp=stdout) const;

        static const double geomRelTol; 
        
        static bool isOrthogonal(const vector3<>&, const vector3<>&);
        
private:
        struct Vertex; 
        struct Edge; 
        struct Face; 

        matrix3<> R, invR, RTR; 
        double minDistSq; 
        std::list<Vertex*> vertex; 
        std::set<Edge*> edge; 
        std::set<Face*> face; 
        std::vector<Face*> faceHalf; 
        
        struct Vertex
        {       vector3<> pos; 
                std::list<Edge*> edge; 
        };
        
        struct Edge
        {       std::array<Vertex*,2> vertex; 
                std::array<Face*,2> face; 
        };
        
        struct Face
        {       vector3<int> img; 
                vector3<> eqn; 
                std::list<Edge*> edge; 
        };
        
        void addPlane(const vector3<int>& a);
        
        void addEdge(Face* f, Vertex* vStart, Vertex* vEnd, bool lastEdge=false);
};

#endif // JDFTX_CORE_WIGNERSEITZ_H