GridInfo.h

Go to the documentation of this file.

/*-------------------------------------------------------------------
Copyright 2011 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_GRIDINFO_H
#define JDFTX_CORE_GRIDINFO_H

#include <core/matrix3.h>
#include <core/ScalarField.h>
#include <core/GpuUtil.h>
#include <map>
#include <fftw3.h>
#include <stdint.h>
#include <cstdio>
#include <mutex>

class GridInfo
{
public:

        GridInfo();
        ~GridInfo();

        void update(); 
        void printLattice(); 
        void printReciprocalLattice(); 
        
        enum LatticeType
        {       Manual, 
                //7 lattice systems specified using sides a,b,c and alpha,beta,gamma
                Triclinic,
                Monoclinic,
                Orthorhombic,
                Tetragonal,
                Rhombohedral,
                Hexagonal,
                Cubic
        }
        latticeType;
        
        enum LatticeModification
        {       Simple,
                BodyCentered,
                BaseCentered,
                FaceCentered
        }
        latticeModification;
        
        double a, b, c, alpha, beta, gamma; 
        void setLatticeVectors(); 
        
        vector3<> lattScale; 
        
        matrix3<> R; 
        double Gmax; 
        double GmaxRho; 
        vector3<int> S; 

        void initialize(bool skipHeader=false, const std::vector< matrix3<int> > sym = std::vector< matrix3<int> >(1, matrix3<int>(1,1,1)));

        double detR; 
        matrix3<> RT, RTR, invR, invRT, invRTR; 
        matrix3<> G, GT, GGT, invGGT; 

        double dV; 
        vector3<> h[3]; 
        int nr; 
        int nG; 

        double dGradial; 
        double GmaxSphere; 
        double GmaxGrid; 
        static const double maxAllowedStrain; 
        
        //Recommended division amongst processes
        //NOTE: Most operators are independent on each process, this is only a suggested division
        //for functions that are expensive to evaluate at each grid point; the corresponding client
        //code is responsible for making sure that the final results are broadcast to all processes
        int irStart, irStop; //division for real space anf full-G space loops
        int iGstart, iGstop; //division for half-G space loops
        
        //FFT plans:
        enum PlanType
        {       PlanForward, 
                PlanInverse, 
                PlanForwardInPlace, 
                PlanInverseInPlace, 
                PlanRtoC, 
                PlanCtoR, 
        };
        fftw_plan getPlan(PlanType planType, int nThreads) const; //get an FFTW plan of specified type with specified thread count
        #ifdef GPU_ENABLED
        cufftHandle planZ2Z; 
        cufftHandle planD2Z; 
        cufftHandle planZ2D; 
        cufftHandle planZ2Dcompat; 
        #endif

        //Indexing utilities (inlined for efficiency)
        inline vector3<int> wrapGcoords(const vector3<int> iG) const 
        {       vector3<int> iGwrapped = iG;
                for(int k=0; k<3; k++)
                        if(iGwrapped[k]<0)
                                iGwrapped[k] += S[k];
                return iGwrapped;
        }
        inline int fullRindex(const vector3<int> iR) const 
        {       return iR[2] + S[2]*(iR[1] + S[1]*iR[0]);
        }
        inline int fullGindex(const vector3<int> iG) const 
        {       return fullRindex(wrapGcoords(iG));
        }
        inline int halfGindex(const vector3<int> iG) const 
        {       vector3<int> iGwrapped = wrapGcoords(iG);
                return iGwrapped[2] + (S[2]/2+1)*(iGwrapped[1] + S[1]*iGwrapped[0]);
        }

private:
        bool initialized; 
        void updateSdependent();
        
        //FFTW plans by thread count and type:
        std::map<std::pair<PlanType,int>,fftw_plan> planCache;
        static std::mutex planLock; //Global lock since planner routines are not thread safe
};

#endif //JDFTX_CORE_GRIDINFO_H