1 /*------------------------------------------------------------------- 2 Copyright 2011 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_GRIDINFO_H 21 #define JDFTX_CORE_GRIDINFO_H 22 23 //! @addtogroup Geometry 24 //! @{ 25 26 #include <core/matrix3.h> 27 #include <core/GpuUtil.h> 28 #include <fftw3.h> 29 #include <stdint.h> 30 #include <cstdio> 31 #include <mutex> 32 #include <map> 33 34 /** @brief Simulation grid descriptor 35 36 To setup a simulation grid, create a blank GridInfo object, 37 set the public data members #S and #R, and call #initialize(). 38 39 This sets up all the auxilliary grid information, and a bunch 40 of shared utilities such as a random number generator, 41 plans for Fourier transforms etc. 42 */ 43 class GridInfo 44 { 45 public: 46 47 GridInfo(); 48 ~GridInfo(); 49 50 void update(); //!< Update the grid information on changing the lattice vectors 51 void printLattice(); //!< Print the lattice vectors 52 void printReciprocalLattice(); //!< Print the reciprocal lattice vectors 53 54 enum LatticeType 55 { Manual, //!< Directly specify R 56 //7 lattice systems specified using sides a,b,c and alpha,beta,gamma 57 Triclinic, 58 Monoclinic, 59 Orthorhombic, 60 Tetragonal, 61 Rhombohedral, 62 Hexagonal, 63 Cubic 64 } 65 latticeType; 66 67 enum LatticeModification 68 { Simple, 69 BodyCentered, 70 BaseCentered, 71 FaceCentered 72 } 73 latticeModification; 74 75 double a, b, c, alpha, beta, gamma; //!< lattice specified by type, modification and standard side lengths, angles in degrees 76 void setLatticeVectors(); //!< Set lattice vectors based on lattice type, modification and parameters above 77 78 vector3<> lattScale; //!< Remember lattice scale specified at input (Note R always includes scale, once latt-scale has processed) 79 80 matrix3<> R; //!< directly specified lattice vectors 81 double Gmax; //!< radius of wavefunction G-sphere, whole density sphere (double the radius) must be inscribable within the FFT box 82 double GmaxRho; //!< if non-zero, override the FFT box inscribable sphere radius 83 vector3<int> S; //!< sample points in each dimension (if 0, will be determined automatically based on Gmax) 84 85 //! Initialize the dependent quantities below. 86 //! If S is specified and is too small for the given Gmax, the call will abort. 87 //! If skipHeader=true, the "initializing the grid" banner will be supressed (useful for auxiliary grids in calculation) 88 //! If sym is specified, then auto-computed S will be ensured commensurate with those symmetries 89 void initialize(bool skipHeader=false, const std::vector<SpaceGroupOp> sym = std::vector<SpaceGroupOp>(1, SpaceGroupOp())); 90 91 double detR; //!< cell volume 92 matrix3<> RT, RTR, invR, invRT, invRTR; //!< various combinations of lattice-vectors 93 matrix3<> G, GT, GGT, invGGT; //!< various combinations of reciporcal lattice-vectors 94 95 double dV; //!< volume per grid point 96 vector3<> h[3]; //!< real space sample vectors 97 int nr; //!< position space grid count = S[0]*S[1]*S[2] 98 int nG; //!< reciprocal lattice count = S[0]*S[1]*(S[2]/2+1) [on account of using r2c and c2r ffts] 99 100 double dGradial; //!< recommended spacing of radial G functions 101 double GmaxSphere; //!< recommended maximum G-vector for radial functions for the wavefunction sphere 102 double GmaxGrid; //!< recommended maximum G-vector for radial functions for the density grid 103 static const double maxAllowedStrain; //!< maximum strain allowed when updating lattice vectors and using update() 104 105 //Recommended division amongst processes 106 //NOTE: Most operators are independent on each process, this is only a suggested division 107 //for functions that are expensive to evaluate at each grid point; the corresponding client 108 //code is responsible for making sure that the final results are broadcast to all processes 109 int irStart, irStop; //division for real space anf full-G space loops 110 int iGstart, iGstop; //division for half-G space loops 111 112 //FFT plans: 113 enum PlanType 114 { PlanForward, //!< Forward complex transform 115 PlanInverse, //!< Inverse complex transform 116 PlanForwardInPlace, //!< Forward in-place complex transform 117 PlanInverseInPlace, //!< Inverse in-place complex transform 118 PlanRtoC, //!< Real to complex transform 119 PlanCtoR, //!< Complex to real transform 120 }; 121 fftw_plan getPlan(PlanType planType, int nThreads) const; //get an FFTW plan of specified type with specified thread count 122 #ifdef GPU_ENABLED 123 cufftHandle planZ2Z; //!< CUFFT plan for all the complex transforms 124 cufftHandle planD2Z; //!< CUFFT plan for R -> G 125 cufftHandle planZ2D; //!< CUFFT plan for G -> R 126 #endif 127 128 //Indexing utilities (inlined for efficiency) wrapGcoords(const vector3<int> iG)129 inline vector3<int> wrapGcoords(const vector3<int> iG) const //!< wrap negative G-indices to the positive side 130 { vector3<int> iGwrapped = iG; 131 for(int k=0; k<3; k++) 132 if(iGwrapped[k]<0) 133 iGwrapped[k] += S[k]; 134 return iGwrapped; 135 } fullRindex(const vector3<int> iR)136 inline int fullRindex(const vector3<int> iR) const //!< Index into the full real-space box: 137 { return iR[2] + S[2]*(iR[1] + S[1]*iR[0]); 138 } fullGindex(const vector3<int> iG)139 inline int fullGindex(const vector3<int> iG) const //!< Index into the full reciprocal-space box: 140 { return fullRindex(wrapGcoords(iG)); 141 } halfGindex(const vector3<int> iG)142 inline int halfGindex(const vector3<int> iG) const //!< Index into the half-reduced reciprocal-space box: 143 { vector3<int> iGwrapped = wrapGcoords(iG); 144 return iGwrapped[2] + (S[2]/2+1)*(iGwrapped[1] + S[1]*iGwrapped[0]); 145 } 146 147 private: 148 bool initialized; //!< keep track of whether initialize() has been called 149 void updateSdependent(); 150 151 //FFTW plans by thread count and type: 152 std::map<std::pair<PlanType,int>,fftw_plan> planCache; 153 static std::mutex planLock; //Global lock since planner routines are not thread safe 154 }; 155 156 //! @} 157 #endif //JDFTX_CORE_GRIDINFO_H 158