LatticeUtils.h

/*-------------------------------------------------------------------
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_LATTICEUTILS_H
#define JDFTX_CORE_LATTICEUTILS_H


#include <core/GridInfo.h>
#include <vector>

static const double symmThreshold = 1e-4;
static const double symmThresholdSq = symmThreshold * symmThreshold;

matrix3<> reduceLatticeVectors(const matrix3<>& R,
        matrix3<int>* transmission=0, matrix3<int>* invTransmission=0);

std::vector<matrix3<int>> getSymmetries(const matrix3<>& R,
        vector3<bool> isTruncated=vector3<bool>(false,false,false),
        matrix3<>* Rreduced=0, matrix3<int>* transmission=0, matrix3<int>* invTransmission=0);

struct Supercell
{
        const GridInfo& gInfo; 
        std::vector<vector3<>> kmesh; 
        matrix3<> Rsuper; 
        matrix3<int> super; 
        
        struct KmeshTransform
        {       unsigned iReduced; 
                unsigned iSym; 
                int invert; 
                vector3<int> offset; 
        };
        std::vector<KmeshTransform> kmeshTransform;
        
        Supercell(const GridInfo& gInfo,
                const std::vector<vector3<>>& kmeshReduced,
                const std::vector<matrix3<int>>& sym, const std::vector<int>& invertList);
};

std::map<vector3<int>, double> getCellMap(const matrix3<>& R, const matrix3<>& Rsup, string fname=string());

//Helper function for PeriodicLookup< vector3<> > used in Supercell::Supercell
inline vector3<> getCoord(const vector3<>& pos) { return pos; }

template<typename T> class PeriodicLookup
{       const std::vector<T>& points;
        vector3<int> S; //lookup mesh sample count
        std::vector< std::vector<size_t> > indices; //list of indices into points, for each  lookup mesh cell
        
        inline size_t meshIndex(vector3<int> iv) const
        {       for(int k=0; k<3; k++) //wrap to [0,S)
                {       if(iv[k] < 0) iv[k] += S[k];
                        if(iv[k] >= S[k]) iv[k] -= S[k];
                }
                return iv[0]+S[0]*size_t(iv[1]+S[1]*iv[2]);
        }

public:
        PeriodicLookup(const std::vector<T>& points, matrix3<> metric, size_t nPointsTarget=0) : points(points)
        {       //Set S such that prod(S) ~ nPoints and the sample counts are proportional to dimension length
                vector3<> Stmp; for(int k=0; k<3; k++) Stmp[k] = sqrt(metric(k,k));
                Stmp *= pow(std::max(points.size(), nPointsTarget)/(Stmp[0]*Stmp[1]*Stmp[2]), 1./3); //normalize so that product is nPoints
                for(int k=0; k<3; k++)
                {       S[k] = std::max(1, int(round(Stmp[k])));
                        assert(symmThreshold*S[k] < 0.5);
                }
                //Initialize indices:
                indices.resize(S[0]*S[1]*S[2]);
                for(size_t iPoint=0; iPoint<points.size(); iPoint++)
                        addPoint(iPoint, points[iPoint]);
        }
        
        void addPoint(size_t iPoint, const T& point)
        {       vector3<> v = getCoord(point);
                vector3<int> iv;
                for(int k=0; k<3; k++)
                {       v[k] -= floor(v[k]); //in [0,1)
                        iv[k] = int(floor(v[k]*S[k] + 0.5)); //in [0, S]
                }
                indices[meshIndex(iv)].push_back(iPoint);
        }
        
        template<typename Tag = double, typename TagEquality = std::equal_to<Tag> >
        size_t find(vector3<> v, Tag tag = Tag(), const std::vector<Tag>* tagArr=0, TagEquality tagEquality = std::equal_to<Tag>()) const
        {       vector3<int> ivMin, ivMax;
                for(int k=0; k<3; k++)
                {       v[k] -= floor(v[k]); //in [0,1)
                        ivMin[k] = int(floor((v[k]-symmThreshold)*S[k] + 0.5)); //in [-1, S]
                        ivMax[k] = int(floor((v[k]+symmThreshold)*S[k] + 0.5)); //in [0, S]
                }
                vector3<int> iv;
                for(iv[0]=ivMin[0]; iv[0]<=ivMax[0]; iv[0]++)
                for(iv[1]=ivMin[1]; iv[1]<=ivMax[1]; iv[1]++)
                for(iv[2]=ivMin[2]; iv[2]<=ivMax[2]; iv[2]++)
                        for(size_t index: indices[meshIndex(iv)])
                                if(circDistanceSquared(v, getCoord(points[index])) < symmThresholdSq
                                        && (!tagArr || tagEquality(tag, tagArr->at(index))) )
                                        return index;
                return string::npos;
        }
};

#endif // JDFTX_CORE_LATTICEUTILS_H