1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 //                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
10 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
11 //                    Ye Luo, yeluo@anl.gov, Argonne National Laboratory
12 //                    Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
13 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
14 //
15 // File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
16 //////////////////////////////////////////////////////////////////////////////////////
17 
18 
19 /** @file EinsplineSetBuilder.h
20  *
21  * Builder class for einspline-based SPOSet objects.
22  */
23 #ifndef QMCPLUSPLUS_EINSPLINE_SET_BUILDER_H
24 #define QMCPLUSPLUS_EINSPLINE_SET_BUILDER_H
25 
26 #include "QMCWaveFunctions/SPOSetBuilder.h"
27 #include "QMCWaveFunctions/BandInfo.h"
28 #include "QMCWaveFunctions/AtomicOrbital.h"
29 #include "QMCWaveFunctions/EinsplineSet.h"
30 #include "Numerics/HDFNumericAttrib.h"
31 #include <map>
32 
33 #define PW_COEFF_NORM_TOLERANCE 1e-6
34 
35 class Communicate;
36 
37 namespace qmcplusplus
38 {
39 ///forward declaration of BsplineReaderBase
40 struct BsplineReaderBase;
41 
42 // Helper needed for TwistMap
43 struct Int3less
44 {
operatorInt3less45   bool operator()(const TinyVector<int, 3>& a, const TinyVector<int, 3>& b) const
46   {
47     if (a[0] > b[0])
48       return false;
49     if (a[0] < b[0])
50       return true;
51     if (a[1] > b[1])
52       return false;
53     if (a[1] < b[1])
54       return true;
55     if (a[2] > b[2])
56       return false;
57     if (a[2] < b[2])
58       return true;
59     return false;
60   }
61 };
62 struct Int4less
63 {
operatorInt4less64   bool operator()(const TinyVector<int, 4>& a, const TinyVector<int, 4>& b) const
65   {
66     for (int i = 0; i < 4; i++)
67     {
68       if (a[i] > b[i])
69         return false;
70       if (a[i] < b[i])
71         return true;
72     }
73     return false;
74   }
75 };
76 
77 
78 /** construct a name for spline SPO set
79  */
80 struct H5OrbSet
81 {
82   ///type of orbitals defined
83   int OrbitalType;
84   ///index for the spin set
85   int SpinSet;
86   ///number of orbitals that belong to this set
87   int NumOrbs;
88   ///name of the HDF5 file
89   std::string FileName;
90   /** true if a < b
91    *
92    * The ordering
93    * - name
94    * - spin set
95    * - number of orbitals
96    */
operatorH5OrbSet97   bool operator()(const H5OrbSet& a, const H5OrbSet& b) const
98   {
99     if (a.FileName == b.FileName)
100     {
101       if (a.SpinSet == b.SpinSet)
102         return a.NumOrbs < b.NumOrbs;
103       else
104         return a.SpinSet < b.SpinSet;
105     }
106     else
107       return a.FileName < b.FileName;
108   }
109 
H5OrbSetH5OrbSet110   H5OrbSet(const H5OrbSet& a) : SpinSet(a.SpinSet), NumOrbs(a.NumOrbs), FileName(a.FileName) {}
H5OrbSetH5OrbSet111   H5OrbSet(std::string name, int spinSet, int numOrbs) : SpinSet(spinSet), NumOrbs(numOrbs), FileName(name) {}
H5OrbSetH5OrbSet112   H5OrbSet() {}
113 };
114 
115 /** EinsplineSet builder
116  */
117 class EinsplineSetBuilder : public SPOSetBuilder
118 {
119 public:
120   typedef std::map<std::string, ParticleSet*> PtclPoolType;
121   typedef CrystalLattice<ParticleSet::Scalar_t, DIM> UnitCellType;
122 
123   ///reference to the particleset pool
124   PtclPoolType& ParticleSets;
125   ///quantum particle set
126   ParticleSet& TargetPtcl;
127   ///ionic system
128   ParticleSet* SourcePtcl;
129 
130   /**  Helper vector for sorting bands
131    */
132   std::vector<std::vector<BandInfo>*> FullBands;
133 
134   /// reader to use BsplineReaderBase
135   BsplineReaderBase* MixedSplineReader;
136 
137   ///This is true if we have the orbital derivatives w.r.t. the ion positions
138   bool HaveOrbDerivs;
139   ///root XML node with href, sort, tilematrix, twistnum, source, precision,truncate,version
140   xmlNodePtr XMLRoot;
141 
142   std::map<H5OrbSet, SPOSet*, H5OrbSet> SPOSetMap;
143 
144   ///constructor
145   EinsplineSetBuilder(ParticleSet& p, PtclPoolType& psets, Communicate* comm, xmlNodePtr cur);
146 
147   ///destructor
148   ~EinsplineSetBuilder();
149 
150   /** initialize the Antisymmetric wave function for electrons
151    * @param cur the current xml node
152    */
153   SPOSet* createSPOSetFromXML(xmlNodePtr cur);
154 
155   /** a specific but clean code path in createSPOSetFromXML, for PBC, double, ESHDF
156    * @param cur the current xml node
157    */
158   void set_metadata(int numOrbs, int TwistNum_inp, bool skipChecks=false);
159 
160   /** initialize with the existing SPOSet */
161   SPOSet* createSPOSet(xmlNodePtr cur, SPOSetInputInfo& input_info);
162 
163   //////////////////////////////////////
164   // HDF5-related data  and functions //
165   //////////////////////////////////////
166   hid_t H5FileID;
167   std::string H5FileName;
168   // HDF5 orbital file version
169   typedef enum
170   {
171     QMCPACK,
172     ESHDF
173   } FormatType;
174   FormatType Format;
175   TinyVector<int, 3> Version;
176   std::string parameterGroup, ionsGroup, eigenstatesGroup;
177   std::vector<int> Occ;
178   bool HasCoreOrbs;
179   bool ReadOrbitalInfo(bool skipChecks=false);
180   bool ReadOrbitalInfo_ESHDF(bool skipChecks=false);
181   void BroadcastOrbitalInfo();
182   bool CheckLattice();
183 
184   /** read gvectors for each twist
185    * @return true, if psi_g is found
186    */
187   bool ReadGvectors_ESHDF();
188 
189   /** set tiling properties of oset
190    * @param oset spline-orbital engine to be initialized
191    * @param numOrbs number of orbitals that belong to oset
192    */
193   template<typename SPE>
setTiling(SPE * oset,int numOrbs)194   inline void setTiling(SPE* oset, int numOrbs)
195   {
196     oset->TileFactor   = TileFactor;
197     oset->Tiling       = (TileFactor[0] * TileFactor[1] * TileFactor[2] != 1);
198     oset->PrimLattice.set(Lattice);
199     oset->SuperLattice.set(SuperLattice);
200     oset->GGt = GGt;
201     oset->setOrbitalSetSize(numOrbs);
202   }
203 
204 
205   Tensor<double, OHMMS_DIM> Lattice, RecipLattice, LatticeInv, SuperLattice, GGt;
206   UnitCellType SuperCell, PrimCell, PrimCellInv;
207   int NumBands, NumElectrons, NumSpins, NumTwists, NumCoreStates;
208   int MaxNumGvecs;
209   double MeshFactor;
210   RealType MatchingTol;
211   TinyVector<int, 3> MeshSize;
212   std::vector<std::vector<TinyVector<int, 3>>> Gvecs;
213 
214   Vector<int> IonTypes;
215   Vector<TinyVector<double, OHMMS_DIM>> IonPos;
216   // mapping the ions in the supercell to the primitive cell
217   std::vector<int> Super2Prim;
218 
219   /////////////////////////////
220   // Twist angle information //
221   /////////////////////////////
222   // This stores which "true" twist number I am using
223   int TwistNum;
224   TinyVector<double, OHMMS_DIM> givenTwist;
225   std::vector<TinyVector<double, OHMMS_DIM>> TwistAngles;
226   //     integer index of sym operation from the irreducible brillion zone
227   std::vector<int> TwistSymmetry;
228   //     number of twists equivalent to this one in the big DFT grid
229   std::vector<int> TwistWeight;
230 
231   TinyVector<int, OHMMS_DIM> TileFactor;
232   Tensor<int, OHMMS_DIM> TileMatrix;
233   TinyVector<int, OHMMS_DIM> TwistMesh;
234   // This vector stores which twist indices will be used by this
235   // clone
236   std::vector<TinyVector<int, OHMMS_DIM>> UseTwists;
237   std::vector<int> IncludeTwists, DistinctTwists;
238   bool UseRealOrbitals;
239   int NumDistinctOrbitals, NumCoreOrbs, NumValenceOrbs;
240   // This is true if the corresponding twist in DistinctTwists should
241   // should be used to generate two distinct orbitals from the real and
242   // imaginary parts.
243   std::vector<bool> MakeTwoCopies;
244   inline bool TwistPair(PosType a, PosType b);
245   // This maps a 3-integer twist index into the twist number in the file
246   std::map<TinyVector<int, OHMMS_DIM>, int, Int3less> TwistMap;
247   //void AnalyzeTwists();
248   void AnalyzeTwists2();
249   void TileIons();
250   void OccupyBands(int spin, int sortBands, int numOrbs, bool skipChecks=false);
251   void OccupyBands_ESHDF(int spin, int sortBands, int numOrbs);
252 
253   void CopyBands(int numOrbs);
254 
255   /////////////////////////////
256   // Muffin-tin information  //
257   /////////////////////////////
258   int NumMuffinTins;
259   std::vector<double> MT_APW_radii;
260   std::vector<Vector<double>> MT_APW_rgrids;
261   std::vector<int> MT_APW_lmax;
262   std::vector<int> MT_APW_num_radial_points;
263   std::vector<TinyVector<double, OHMMS_DIM>> MT_centers;
264 
265   ////////////////////////////////
266   // Atomic orbital information //
267   ////////////////////////////////
268   std::vector<AtomicOrbital<std::complex<double>>> AtomicOrbitals;
269 
270   struct CenterInfo
271   {
272     std::vector<int> lmax, spline_npoints, GroupID;
273     std::vector<double> spline_radius, cutoff, inner_cutoff, non_overlapping_radius;
274     std::vector<TinyVector<double, OHMMS_DIM>> ion_pos;
275     int Ncenters;
276 
CenterInfoCenterInfo277     CenterInfo() : Ncenters(0){};
278 
resizeCenterInfo279     void resize(int ncenters)
280     {
281       Ncenters = ncenters;
282       lmax.resize(ncenters, -1);
283       spline_npoints.resize(ncenters, -1);
284       GroupID.resize(ncenters, 0);
285       spline_radius.resize(ncenters, -1.0);
286       inner_cutoff.resize(ncenters, -1.0);
287       non_overlapping_radius.resize(ncenters, -1.0);
288       cutoff.resize(ncenters, -1.0);
289       ion_pos.resize(ncenters);
290     }
291   } AtomicCentersInfo;
292 
293   // This returns the path in the HDF5 file to the group for orbital
294   // with twist ti and band bi
295   std::string OrbitalPath(int ti, int bi);
296   std::string CoreStatePath(int ti, int bi);
297   std::string MuffinTinPath(int ti, int bi, int tin);
298 
299   /////////////////////////////////////////////////////////////
300   // Information to avoid storing the same orbitals twice in //
301   // spin-restricted calculations.                           //
302   /////////////////////////////////////////////////////////////
303   int LastSpinSet, NumOrbitalsRead;
304 
305   std::string occ_format;
306   int particle_hole_pairs;
307   bool makeRotations;
308 
309   /** broadcast SortBands
310    * @param N number of state
311    * @param root true if it is the i/o node
312    * @return true, if core is found
313    */
314   bool bcastSortBands(int splin, int N, bool root);
315 
316   int MyToken;
update_token(const char * f,int l,const char * msg)317   inline void update_token(const char* f, int l, const char* msg)
318   {
319     app_debug() << "TOKEN=" << MyToken << " " << msg << " " << f << " " << l << std::endl;
320     MyToken++;
321   }
322   //inline void update_token(const char* f, int l, const char* msg)
323   //{}
324 };
325 
326 } // namespace qmcplusplus
327 
328 
329 #endif
330