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 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 // Ye Luo, yeluo@anl.gov, Argonne National Laboratory
11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //
13 // File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15
16
17 #include "QMCWaveFunctions/EinsplineSetBuilder.h"
18 #include "QMCWaveFunctions/WaveFunctionComponentBuilder.h"
19 #include "OhmmsData/AttributeSet.h"
20 #include "Utilities/Timer.h"
21 #include "Message/Communicate.h"
22 #include "Message/CommOperators.h"
23 #include <vector>
24 #include "Numerics/HDFSTLAttrib.h"
25 #include "OhmmsData/HDFStringAttrib.h"
26 #include "ParticleIO/ESHDFParticleParser.h"
27 #include "ParticleBase/RandomSeqGenerator.h"
28 #include "config/stdlib/math.hpp"
29
30 namespace qmcplusplus
31 {
ReadOrbitalInfo(bool skipChecks)32 bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks)
33 {
34 update_token(__FILE__, __LINE__, "ReadOrbitalInfo");
35 // Handle failed file open gracefully by temporarily replacing error handler
36 H5E_auto2_t old_efunc;
37 void* old_efunc_data;
38 H5Eget_auto2(H5E_DEFAULT, &old_efunc, &old_efunc_data);
39 H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
40 H5FileID = H5Fopen(H5FileName.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
41 // H5FileID = H5Fopen(H5FileName.c_str(),H5F_ACC_RDWR,H5P_DEFAULT);
42 if (H5FileID < 0)
43 {
44 app_error() << "Could not open HDF5 file \"" << H5FileName
45 << "\" in EinsplineSetBuilder::ReadOrbitalInfo. Aborting.\n";
46 APP_ABORT("EinsplineSetBuilder::ReadOrbitalInfo");
47 }
48 H5Eset_auto2(H5E_DEFAULT, old_efunc, old_efunc_data);
49
50 // Read format
51 std::string format;
52 HDFAttribIO<std::string> h_format(format);
53 h_format.read(H5FileID, "/format");
54 HDFAttribIO<TinyVector<int, 3>> h_Version(Version);
55 h_Version.read(H5FileID, "/version");
56 app_log() << " HDF5 orbital file version " << Version[0] << "." << Version[1] << "." << Version[2] << "\n";
57 if (format.find("ES") < format.size())
58 {
59 Format = ESHDF;
60 return ReadOrbitalInfo_ESHDF(skipChecks);
61 }
62 //////////////////////////////////////////////////
63 // Read basic parameters from the orbital file. //
64 //////////////////////////////////////////////////
65 // Check the version
66 if (Version[0] == 0 && Version[1] == 11)
67 {
68 parameterGroup = "/parameters_0";
69 ionsGroup = "/ions_2";
70 eigenstatesGroup = "/eigenstates_3";
71 }
72 else if (Version[0] == 0 && Version[1] == 20)
73 {
74 parameterGroup = "/parameters";
75 ionsGroup = "/ions";
76 eigenstatesGroup = "/eigenstates";
77 }
78 else
79 {
80 std::ostringstream o;
81 o << "Unknown HDF5 orbital file version " << Version[0] << "." << Version[1] << "." << Version[2] << "\n";
82 APP_ABORT(o.str());
83 //abort();
84 }
85 HDFAttribIO<Tensor<double, 3>> h_Lattice(Lattice), h_RecipLattice(RecipLattice);
86 h_Lattice.read(H5FileID, (parameterGroup + "/lattice").c_str());
87 h_RecipLattice.read(H5FileID, (parameterGroup + "/reciprocal_lattice").c_str());
88 SuperLattice = dot(TileMatrix, Lattice);
89 char buff[1000];
90 snprintf(buff, 1000,
91 " Lattice = \n [ %8.5f %8.5f %8.5f\n"
92 " %8.5f %8.5f %8.5f\n"
93 " %8.5f %8.5f %8.5f ]\n",
94 Lattice(0, 0), Lattice(0, 1), Lattice(0, 2), Lattice(1, 0), Lattice(1, 1), Lattice(1, 2), Lattice(2, 0),
95 Lattice(2, 1), Lattice(2, 2));
96 app_log() << buff;
97 snprintf(buff, 1000,
98 " SuperLattice = \n [ %13.12f %13.12f %13.12f\n"
99 " %13.12f %13.12f %13.12f\n"
100 " %13.12f %13.12f %13.12f ]\n",
101 SuperLattice(0, 0), SuperLattice(0, 1), SuperLattice(0, 2), SuperLattice(1, 0), SuperLattice(1, 1),
102 SuperLattice(1, 2), SuperLattice(2, 0), SuperLattice(2, 1), SuperLattice(2, 2));
103 if (!CheckLattice())
104 APP_ABORT("CheckLattice failed");
105 app_log() << buff;
106 for (int i = 0; i < 3; i++)
107 for (int j = 0; j < 3; j++)
108 LatticeInv(i, j) = RecipLattice(i, j) / (2.0 * M_PI);
109 HDFAttribIO<int> h_NumBands(NumBands), h_NumElectrons(NumElectrons), h_NumSpins(NumSpins), h_NumTwists(NumTwists),
110 h_NumCore(NumCoreStates), h_NumMuffinTins(NumMuffinTins);
111 NumCoreStates = NumMuffinTins = 0;
112 h_NumBands.read(H5FileID, (parameterGroup + "/num_bands").c_str());
113 h_NumCore.read(H5FileID, (parameterGroup + "/num_core_states").c_str());
114 h_NumElectrons.read(H5FileID, (parameterGroup + "/num_electrons").c_str());
115 h_NumSpins.read(H5FileID, (parameterGroup + "/num_spins").c_str());
116 h_NumTwists.read(H5FileID, (parameterGroup + "/num_twists").c_str());
117 h_NumMuffinTins.read(H5FileID, (parameterGroup + "/muffin_tins/num_tins").c_str());
118 app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists
119 << ", muffin tins=" << NumMuffinTins << std::endl;
120 // fprintf (stderr, " bands = %d, elecs = %d, spins = %d, twists = %d\n",
121 // NumBands, NumElectrons, NumSpins, NumTwists);
122 if (TileFactor[0] != 1 || TileFactor[1] != 1 || TileFactor[2] != 1)
123 app_log() << " Using a " << TileFactor[0] << "x" << TileFactor[1] << "x" << TileFactor[2] << " tiling factor.\n";
124 /////////////////////////////////
125 // Read muffin tin information //
126 /////////////////////////////////
127 MT_APW_radii.resize(NumMuffinTins);
128 MT_APW_rgrids.resize(NumMuffinTins);
129 MT_APW_lmax.resize(NumMuffinTins);
130 MT_APW_num_radial_points.resize(NumMuffinTins);
131 MT_centers.resize(NumMuffinTins);
132 for (int tin = 0; tin < NumMuffinTins; tin++)
133 {
134 std::ostringstream MTstream;
135 if (NumMuffinTins > 1)
136 MTstream << parameterGroup << "/muffin_tins/muffin_tin_" << tin;
137 else
138 MTstream << parameterGroup << "/muffin_tins/muffin_tin";
139 std::string MTgroup = MTstream.str();
140 HDFAttribIO<int> h_lmax(MT_APW_lmax[tin]), h_num_radial_points(MT_APW_num_radial_points[tin]);
141 HDFAttribIO<double> h_radius(MT_APW_radii[tin]);
142 HDFAttribIO<TinyVector<double, OHMMS_DIM>> h_center(MT_centers[tin]);
143 HDFAttribIO<Vector<double>> h_rgrid(MT_APW_rgrids[tin]);
144 h_lmax.read(H5FileID, (MTgroup + "/lmax").c_str());
145 h_num_radial_points.read(H5FileID, (MTgroup + "/num_radial_points").c_str());
146 h_radius.read(H5FileID, (MTgroup + "/radius").c_str());
147 h_center.read(H5FileID, (MTgroup + "/center").c_str());
148 h_rgrid.read(H5FileID, (MTgroup + "/r").c_str());
149 }
150 //////////////////////////////////
151 // Read ion types and locations //
152 //////////////////////////////////
153 HDFAttribIO<Vector<int>> h_IonTypes(IonTypes);
154 HDFAttribIO<Vector<TinyVector<double, 3>>> h_IonPos(IonPos);
155 h_IonTypes.read(H5FileID, (ionsGroup + "/atom_types").c_str());
156 h_IonPos.read(H5FileID, (ionsGroup + "/pos").c_str());
157 ///////////////////////////
158 // Read the twist angles //
159 ///////////////////////////
160 TwistAngles.resize(NumTwists);
161 for (int ti = 0; ti < NumTwists; ti++)
162 {
163 std::ostringstream path;
164 if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1)
165 path << eigenstatesGroup << "/twist_" << ti << "/twist_angle";
166 else
167 path << eigenstatesGroup << "/twist/twist_angle";
168 TinyVector<double, OHMMS_DIM> TwistAngles_DP;
169 HDFAttribIO<TinyVector<double, OHMMS_DIM>> h_Twist(TwistAngles_DP);
170 h_Twist.read(H5FileID, path.str().c_str());
171 TwistAngles[ti] = TwistAngles_DP;
172 snprintf(buff, 1000, " Found twist angle (%6.3f, %6.3f, %6.3f)\n", TwistAngles[ti][0], TwistAngles[ti][1],
173 TwistAngles[ti][2]);
174 app_log() << buff;
175 }
176 //////////////////////////////////////////////////////////
177 // If the density has not been set in TargetPtcl, and //
178 // the density is available, read it in and save it //
179 // in TargetPtcl. //
180 //////////////////////////////////////////////////////////
181 if (!TargetPtcl.Density_G.size())
182 {
183 HDFAttribIO<std::vector<TinyVector<int, OHMMS_DIM>>> h_reduced_gvecs(TargetPtcl.DensityReducedGvecs);
184 Array<double, OHMMS_DIM> Density_r_DP;
185 HDFAttribIO<Array<double, OHMMS_DIM>> h_density_r(Density_r_DP);
186 h_reduced_gvecs.read(H5FileID, "/density/reduced_gvecs");
187 h_density_r.read(H5FileID, "/density/rho_r");
188 TargetPtcl.Density_r = Density_r_DP;
189 int numG = TargetPtcl.DensityReducedGvecs.size();
190 // Convert primitive G-vectors to supercell G-vectors
191 for (int iG = 0; iG < numG; iG++)
192 TargetPtcl.DensityReducedGvecs[iG] = dot(TileMatrix, TargetPtcl.DensityReducedGvecs[iG]);
193 app_log() << " Read " << numG << " density G-vectors.\n";
194 if (TargetPtcl.DensityReducedGvecs.size())
195 {
196 app_log() << " EinsplineSetBuilder found density in the HDF5 file.\n";
197 std::vector<std::complex<double>> Density_G_DP;
198 HDFAttribIO<std::vector<std::complex<double>>> h_density_G(Density_G_DP);
199 h_density_G.read(H5FileID, "/density/rho_G");
200 TargetPtcl.Density_G.assign(Density_G_DP.begin(), Density_G_DP.end());
201 if (!TargetPtcl.Density_G.size())
202 {
203 app_error() << " Density reduced G-vectors defined, but not the"
204 << " density.\n";
205 abort();
206 }
207 }
208 }
209 return true;
210 }
211
212
OrbitalPath(int ti,int bi)213 std::string EinsplineSetBuilder::OrbitalPath(int ti, int bi)
214 {
215 std::string eigenstatesGroup;
216 if (Version[0] == 0 && Version[1] == 11)
217 eigenstatesGroup = "/eigenstates_3";
218 else if (Version[0] == 0 && Version[1] == 20)
219 eigenstatesGroup = "/eigenstates";
220 std::ostringstream groupPath;
221 if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1)
222 groupPath << eigenstatesGroup << "/twist_" << ti << "/band_" << bi << "/";
223 else if (NumBands > 1)
224 groupPath << eigenstatesGroup << "/twist/band_" << bi << "/";
225 else
226 groupPath << eigenstatesGroup << "/twist/band/";
227 return groupPath.str();
228 }
229
CoreStatePath(int ti,int cs)230 std::string EinsplineSetBuilder::CoreStatePath(int ti, int cs)
231 {
232 std::string eigenstatesGroup;
233 if (Version[0] == 0 && Version[1] == 11)
234 eigenstatesGroup = "/eigenstates_3";
235 else if (Version[0] == 0 && Version[1] == 20)
236 eigenstatesGroup = "/eigenstates";
237 std::ostringstream groupPath;
238 if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1)
239 groupPath << eigenstatesGroup << "/twist_" << ti << "/core_state_" << cs << "/";
240 else if (NumBands > 1)
241 groupPath << eigenstatesGroup << "/twist/core_state_" << cs << "/";
242 else
243 groupPath << eigenstatesGroup << "/twist/core_state/";
244 return groupPath.str();
245 }
246
MuffinTinPath(int ti,int bi,int tin)247 std::string EinsplineSetBuilder::MuffinTinPath(int ti, int bi, int tin)
248 {
249 std::ostringstream groupPath;
250 if (NumMuffinTins > 0)
251 groupPath << OrbitalPath(ti, bi) << "muffin_tin_" << tin << "/";
252 else
253 groupPath << OrbitalPath(ti, bi) << "muffin_tin/";
254 return groupPath.str();
255 }
256
257 #if 0
258 void
259 EinsplineSetBuilder::ReadBands
260 (int spin, EinsplineSetExtended<std::complex<double> >* orbitalSet)
261 {
262 update_token(__FILE__,__LINE__,"ReadBands:complex");
263 bool root = myComm->rank()==0;
264 //bcastwith other stuff
265 myComm->bcast(NumDistinctOrbitals);
266 myComm->bcast (NumValenceOrbs);
267 myComm->bcast (NumCoreOrbs);
268 int N = NumDistinctOrbitals;
269 orbitalSet->kPoints.resize(N);
270 orbitalSet->MakeTwoCopies.resize(N);
271 orbitalSet->StorageValueVector.resize(N);
272 orbitalSet->BlendValueVector.resize(N);
273 orbitalSet->StorageLaplVector.resize(N);
274 orbitalSet->BlendLaplVector.resize(N);
275 orbitalSet->StorageGradVector.resize(N);
276 orbitalSet->BlendGradVector.resize(N);
277 orbitalSet->StorageHessVector.resize(N);
278 orbitalSet->phase.resize(N);
279 orbitalSet->eikr.resize(N);
280 orbitalSet->NumValenceOrbs = NumValenceOrbs;
281 orbitalSet->NumCoreOrbs = NumCoreOrbs;
282
283 std::vector<BandInfo>& SortBands(*FullBands[spin]);
284 if (root)
285 {
286 int numOrbs = orbitalSet->getOrbitalSetSize();
287 int num = 0;
288 for (int iorb=0; iorb<N; iorb++)
289 {
290 int ti = SortBands[iorb].TwistIndex;
291 PosType twist = TwistAngles[ti];
292 orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
293 orbitalSet->MakeTwoCopies[iorb] =
294 (num < (numOrbs-1)) && SortBands[iorb].MakeTwoCopies;
295 num += orbitalSet->MakeTwoCopies[iorb] ? 2 : 1;
296 }
297 }
298 myComm->bcast(orbitalSet->kPoints);
299 myComm->bcast(orbitalSet->MakeTwoCopies);
300 // First, check to see if we have already read this in
301 H5OrbSet set(H5FileName, spin, N);
302 // std::map<H5OrbSet,multi_UBspline_3d_z*>::iterator iter;
303 // iter = ExtendedMap_z.find (set);
304 // if (iter != ExtendedMap_z.end()) {
305 // std::cerr << "Using existing copy of multi_UBspline_3d_z for "
306 // << "thread number " << omp_get_thread_num() << ".\n";
307 // orbitalSet->MultiSpline = iter->second;
308 // return;
309 // }
310 int nx, ny, nz, bi, ti;
311 Array<std::complex<double>,3> splineData, rawData;
312 if (root)
313 {
314 // Find the orbital mesh size
315 int i=0;
316 while (SortBands[i].IsCoreState)
317 i++;
318 ti = SortBands[i].TwistIndex;
319 bi = SortBands[i].BandIndex;
320 std::string vectorName = OrbitalPath (ti, bi) + "eigenvector";
321 HDFAttribIO<Array<std::complex<double>,3> > h_rawData(rawData);
322 h_rawData.read(H5FileID, vectorName.c_str());
323 nx = rawData.size(0);
324 ny = rawData.size(1);
325 nz = rawData.size(2);
326 splineData.resize(nx-1, ny-1, nz-1);
327 for (int ix=0; ix<(nx-1); ix++)
328 for (int iy=0; iy<(ny-1); iy++)
329 for (int iz=0; iz<(nz-1); iz++)
330 splineData(ix,iy,iz) = rawData(ix,iy,iz);
331 PosType twist, k;
332 twist = TwistAngles[ti];
333 k = orbitalSet->PrimLattice.k_cart(twist);
334 double e = SortBands[i].Energy;
335 // fprintf (stderr, " ti=%3d bi=%3d energy=%8.5f k=(%7.4f, %7.4f, %7.4f) rank=%d\n",
336 // ti, bi, e, k[0], k[1], k[2], myComm->rank());
337 }
338 TinyVector<int,3> nxyz(nx,ny,nz);
339 myComm->bcast(nxyz);
340 nx=nxyz[0];
341 ny=nxyz[1];
342 nz=nxyz[2];
343 if (!root)
344 splineData.resize(nx-1,ny-1,nz-1);
345 myComm->bcast(splineData);
346 Ugrid x_grid, y_grid, z_grid;
347 BCtype_z xBC, yBC, zBC;
348 xBC.lCode = PERIODIC;
349 xBC.rCode = PERIODIC;
350 yBC.lCode = PERIODIC;
351 yBC.rCode = PERIODIC;
352 zBC.lCode = PERIODIC;
353 zBC.rCode = PERIODIC;
354 x_grid.start = 0.0;
355 x_grid.end = 1.0;
356 x_grid.num = nx-1;
357 y_grid.start = 0.0;
358 y_grid.end = 1.0;
359 y_grid.num = ny-1;
360 z_grid.start = 0.0;
361 z_grid.end = 1.0;
362 z_grid.num = nz-1;
363 // Create the multiUBspline object
364 orbitalSet->MultiSpline =
365 create_multi_UBspline_3d_z (x_grid, y_grid, z_grid, xBC, yBC, zBC, NumValenceOrbs);
366 set_multi_UBspline_3d_z (orbitalSet->MultiSpline, 0, splineData.data());
367 //////////////////////////////////////
368 // Create the MuffinTin APW splines //
369 //////////////////////////////////////
370 orbitalSet->MuffinTins.resize(NumMuffinTins);
371 for (int tin=0; tin<NumMuffinTins; tin++)
372 {
373 orbitalSet->MuffinTins[tin].Atom = tin;
374 orbitalSet->MuffinTins[tin].set_center (MT_centers[tin]);
375 orbitalSet->MuffinTins[tin].set_lattice(Lattice);
376 orbitalSet->MuffinTins[tin].init_APW
377 (MT_APW_rgrids[tin], MT_APW_lmax[tin],
378 NumValenceOrbs);
379 }
380 int iorb = 0;
381 int icore = 0;
382 int ival = 0;
383 while (iorb < N)
384 {
385 bool isCore;
386 if (root)
387 isCore = SortBands[iorb].IsCoreState;
388 myComm->bcast (isCore);
389 if (isCore)
390 {
391 int atom, l, m=0;
392 double rmax;
393 Vector<double> g, r;
394 PosType twist, k;
395 if (root)
396 {
397 ti = SortBands[iorb].TwistIndex;
398 bi = SortBands[iorb].BandIndex;
399 double e = SortBands[iorb].Energy;
400 twist = TwistAngles[ti];
401 k = orbitalSet->PrimLattice.k_cart(twist);
402 std::string atomName = CoreStatePath (ti, bi) + "atom";
403 std::string gName = CoreStatePath (ti, bi) + "g";
404 std::string rMaxName = CoreStatePath (ti, bi) + "rmax";
405 std::string lName = CoreStatePath (ti, bi) + "l";
406 std::string kName = CoreStatePath (ti, bi) + "k";
407 std::string rName = CoreStatePath (ti, bi) + "r";
408 HDFAttribIO<int> h_atom(atom), h_l(l);
409 HDFAttribIO<double> h_rmax(rmax);
410 HDFAttribIO<Vector<double> > h_g(g);
411 HDFAttribIO<Vector<double> > h_r(r);
412 h_atom.read (H5FileID, atomName.c_str());
413 h_l.read (H5FileID, lName.c_str());
414 h_rmax.read (H5FileID, rMaxName.c_str());
415 h_g.read (H5FileID, gName.c_str());
416 h_r.read (H5FileID, rName.c_str());
417 fprintf (stderr, " Core state: ti=%3d bi=%3d energy=%8.5f "
418 "k=(%7.4f, %7.4f, %7.4f) rank=%d\n",
419 ti, bi, e, k[0], k[1], k[2], myComm->rank());
420 }
421 myComm->bcast (atom);
422 myComm->bcast(rmax);
423 myComm->bcast (l);
424 myComm->bcast (k);
425 int ng = g.size();
426 myComm->bcast(ng);
427 if (g.size() != ng)
428 {
429 g.resize(ng);
430 r.resize(ng);
431 }
432 myComm->bcast (g);
433 myComm->bcast (r);
434 double Z = (double)IonTypes[atom];
435 orbitalSet->MuffinTins[atom].addCore (l, m, r, g, k, Z);
436 icore++;
437 }
438 else
439 {
440 if (root)
441 {
442 int ti = SortBands[iorb].TwistIndex;
443 int bi = SortBands[iorb].BandIndex;
444 double e = SortBands[iorb].Energy;
445 PosType twist, k;
446 twist = TwistAngles[ti];
447 k = orbitalSet->PrimLattice.k_cart(twist);
448 char vs[256];
449 sprintf (vs, " Valence state: ti=%3d bi=%3d energy=%8.5f k=(%7.4f, %7.4f, %7.4f) rank=%d\n",
450 ti, bi, e, k[0], k[1], k[2], myComm->rank());
451 app_log() << vs << std::endl;
452
453 std::string vectorName = OrbitalPath (ti, bi) + "eigenvector";
454 HDFAttribIO<Array<std::complex<double>,3> > h_rawData(rawData);
455 h_rawData.read(H5FileID, vectorName.c_str());
456 if ((rawData.size(0) != nx) ||
457 (rawData.size(1) != ny) ||
458 (rawData.size(2) != nz))
459 {
460 fprintf (stderr, "Error in EinsplineSetBuilder::ReadBands.\n");
461 fprintf (stderr, "Extended orbitals should all have the same dimensions\n");
462 abort();
463 }
464 for (int ix=0; ix<(nx-1); ix++)
465 for (int iy=0; iy<(ny-1); iy++)
466 for (int iz=0; iz<(nz-1); iz++)
467 splineData(ix,iy,iz) = rawData(ix,iy,iz);
468 }
469 myComm->bcast(splineData);
470 set_multi_UBspline_3d_z
471 (orbitalSet->MultiSpline, ival, splineData.data());
472 // Now read muffin tin data
473 for (int tin=0; tin<NumMuffinTins; tin++)
474 {
475 // app_log() << "Reading data for muffin tin " << tin << std::endl;
476 PosType twist, k;
477 int lmax = MT_APW_lmax[tin];
478 int numYlm = (lmax+1)*(lmax+1);
479 Array<std::complex<double>,2>
480 u_lm_r(numYlm, MT_APW_num_radial_points[tin]);
481 Array<std::complex<double>,1> du_lm_dr (numYlm);
482 if (root)
483 {
484 int ti = SortBands[iorb].TwistIndex;
485 int bi = SortBands[iorb].BandIndex;
486 twist = TwistAngles[ti];
487 k = orbitalSet->PrimLattice.k_cart(twist);
488 std::string uName = MuffinTinPath (ti, bi,tin) + "u_lm_r";
489 std::string duName = MuffinTinPath (ti, bi,tin) + "du_lm_dr";
490 HDFAttribIO<Array<std::complex<double>,2> > h_u_lm_r(u_lm_r);
491 HDFAttribIO<Array<std::complex<double>,1> > h_du_lm_dr(du_lm_dr);
492 h_u_lm_r.read(H5FileID, uName.c_str());
493 h_du_lm_dr.read(H5FileID, duName.c_str());
494 }
495 myComm->bcast(u_lm_r);
496 myComm->bcast(du_lm_dr);
497 myComm->bcast(k);
498 double Z = (double)IonTypes[tin];
499 orbitalSet->MuffinTins[tin].set_APW (ival, k, u_lm_r, du_lm_dr, Z);
500 }
501 ival++;
502 } // valence state
503 iorb++;
504 }
505 //ExtendedMap_z[set] = orbitalSet->MultiSpline;
506 }
507
508 void
509 EinsplineSetBuilder::ReadBands
510 (int spin, EinsplineSetExtended<double>* orbitalSet)
511 {
512 update_token(__FILE__,__LINE__,"ReadBands:double");
513 std::vector<BandInfo>& SortBands(*FullBands[spin]);
514 bool root = myComm->rank()==0;
515 // bcast other stuff
516 myComm->bcast (NumDistinctOrbitals);
517 myComm->bcast (NumValenceOrbs);
518 myComm->bcast (NumCoreOrbs);
519 int N = NumDistinctOrbitals;
520 orbitalSet->kPoints.resize(N);
521 orbitalSet->MakeTwoCopies.resize(N);
522 orbitalSet->StorageValueVector.resize(N);
523 orbitalSet->BlendValueVector.resize(N);
524 orbitalSet->StorageLaplVector.resize(N);
525 orbitalSet->BlendLaplVector.resize(N);
526 orbitalSet->StorageGradVector.resize(N);
527 orbitalSet->BlendGradVector.resize(N);
528 orbitalSet->StorageHessVector.resize(N);
529 orbitalSet->phase.resize(N);
530 orbitalSet->eikr.resize(N);
531 orbitalSet->NumValenceOrbs = NumValenceOrbs;
532 orbitalSet->NumCoreOrbs = NumCoreOrbs;
533 // Read in k-points
534 int numOrbs = orbitalSet->getOrbitalSetSize();
535 int num = 0;
536 if (root)
537 {
538 for (int iorb=0; iorb<N; iorb++)
539 {
540 int ti = SortBands[iorb].TwistIndex;
541 PosType twist = TwistAngles[ti];
542 orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
543 orbitalSet->MakeTwoCopies[iorb] =
544 (num < (numOrbs-1)) && SortBands[iorb].MakeTwoCopies;
545 num += orbitalSet->MakeTwoCopies[iorb] ? 2 : 1;
546 }
547 PosType twist0 = TwistAngles[SortBands[0].TwistIndex];
548 for (int i=0; i<OHMMS_DIM; i++)
549 if (std::abs(std::abs(twist0[i]) - 0.5) < 1.0e-8)
550 orbitalSet->HalfG[i] = 1;
551 else
552 orbitalSet->HalfG[i] = 0;
553 }
554 myComm->bcast(orbitalSet->kPoints);
555 myComm->bcast(orbitalSet->MakeTwoCopies);
556 myComm->bcast(orbitalSet->HalfG);
557 // First, check to see if we have already read this in
558 H5OrbSet set(H5FileName, spin, N);
559 int nx, ny, nz, bi, ti;
560 Array<std::complex<double>,3> rawData;
561 Array<double,3> splineData;
562 if (root)
563 {
564 // Find the orbital mesh size
565 int i=0;
566 while (SortBands[i].IsCoreState)
567 i++;
568 ti = SortBands[i].TwistIndex;
569 bi = SortBands[i].BandIndex;
570 std::string vectorName = OrbitalPath (ti, bi) + "eigenvector";
571 HDFAttribIO<Array<std::complex<double>,3> > h_rawData(rawData);
572 h_rawData.read(H5FileID, vectorName.c_str());
573 nx = rawData.size(0);
574 ny = rawData.size(1);
575 nz = rawData.size(2);
576 splineData.resize(nx-1, ny-1, nz-1);
577 PosType ru;
578 for (int ix=0; ix<(nx-1); ix++)
579 {
580 ru[0] = (RealType)ix / (RealType)(nx-1);
581 for (int iy=0; iy<(ny-1); iy++)
582 {
583 ru[1] = (RealType)iy / (RealType)(ny-1);
584 for (int iz=0; iz<(nz-1); iz++)
585 {
586 ru[2] = (RealType)iz / (RealType)(nz-1);
587 double phi = -2.0*M_PI*dot (ru, TwistAngles[ti]);
588 double s, c;
589 qmcplusplus::sincos(phi, &s, &c);
590 std::complex<double> phase(c,s);
591 std::complex<double> z = phase*rawData(ix,iy,iz);
592 splineData(ix,iy,iz) = z.imag();
593 }
594 }
595 }
596 PosType twist, k;
597 twist = TwistAngles[ti];
598 k = orbitalSet->PrimLattice.k_cart(twist);
599 double e = SortBands[i].Energy;
600 }
601 TinyVector<int,3> nxyz(nx,ny,nz);
602 myComm->bcast(nxyz);
603 nx=nxyz[0];
604 ny=nxyz[1];
605 nz=nxyz[2];
606 if (!root)
607 splineData.resize(nx-1,ny-1,nz-1);
608 myComm->bcast(splineData);
609 Ugrid x_grid, y_grid, z_grid;
610 BCtype_d xBC, yBC, zBC;
611 if (orbitalSet->HalfG[0])
612 {
613 xBC.lCode = ANTIPERIODIC;
614 xBC.rCode = ANTIPERIODIC;
615 }
616 else
617 {
618 xBC.lCode = PERIODIC;
619 xBC.rCode = PERIODIC;
620 }
621 if (orbitalSet->HalfG[1])
622 {
623 yBC.lCode = ANTIPERIODIC;
624 yBC.rCode = ANTIPERIODIC;
625 }
626 else
627 {
628 yBC.lCode = PERIODIC;
629 yBC.rCode = PERIODIC;
630 }
631 if (orbitalSet->HalfG[2])
632 {
633 zBC.lCode = ANTIPERIODIC;
634 zBC.rCode = ANTIPERIODIC;
635 }
636 else
637 {
638 zBC.lCode = PERIODIC;
639 zBC.rCode = PERIODIC;
640 }
641 x_grid.start = 0.0;
642 x_grid.end = 1.0;
643 x_grid.num = nx-1;
644 y_grid.start = 0.0;
645 y_grid.end = 1.0;
646 y_grid.num = ny-1;
647 z_grid.start = 0.0;
648 z_grid.end = 1.0;
649 z_grid.num = nz-1;
650 // Create the multiUBspline object
651 orbitalSet->MultiSpline =
652 create_multi_UBspline_3d_d (x_grid, y_grid, z_grid, xBC, yBC, zBC, NumValenceOrbs);
653 set_multi_UBspline_3d_d (orbitalSet->MultiSpline, 0, splineData.data());
654 //////////////////////////////////////
655 // Create the MuffinTin APW splines //
656 //////////////////////////////////////
657 orbitalSet->MuffinTins.resize(NumMuffinTins);
658 for (int tin=0; tin<NumMuffinTins; tin++)
659 {
660 orbitalSet->MuffinTins[tin].Atom = tin;
661 orbitalSet->MuffinTins[tin].set_center (MT_centers[tin]);
662 orbitalSet->MuffinTins[tin].set_lattice(Lattice);
663 orbitalSet->MuffinTins[tin].init_APW
664 (MT_APW_rgrids[tin], MT_APW_lmax[tin],
665 NumValenceOrbs);
666 }
667 int iorb = 0;
668 int icore = 0;
669 int ival = 0;
670 while (iorb < N)
671 {
672 bool isCore;
673 if (root)
674 isCore = SortBands[iorb].IsCoreState;
675 myComm->bcast (isCore);
676 if (isCore)
677 {
678 int atom, l, m=0;
679 double rmax;
680 Vector<double> g, r;
681 PosType twist, k;
682 if (root)
683 {
684 ti = SortBands[iorb].TwistIndex;
685 bi = SortBands[iorb].BandIndex;
686 double e = SortBands[iorb].Energy;
687 twist = TwistAngles[ti];
688 k = orbitalSet->PrimLattice.k_cart(twist);
689 std::string atomName = CoreStatePath (ti, bi) + "atom";
690 std::string gName = CoreStatePath (ti, bi) + "g";
691 std::string rMaxName = CoreStatePath (ti, bi) + "rmax";
692 std::string lName = CoreStatePath (ti, bi) + "l";
693 std::string kName = CoreStatePath (ti, bi) + "k";
694 std::string rName = CoreStatePath (ti, bi) + "r";
695 HDFAttribIO<int> h_atom(atom), h_l(l);
696 HDFAttribIO<double> h_rmax(rmax);
697 HDFAttribIO<Vector<double> > h_g(g);
698 HDFAttribIO<Vector<double> > h_r(r);
699 h_atom.read (H5FileID, atomName.c_str());
700 h_l.read (H5FileID, lName.c_str());
701 h_rmax.read (H5FileID, rMaxName.c_str());
702 h_g.read (H5FileID, gName.c_str());
703 h_r.read (H5FileID, rName.c_str());
704 fprintf (stderr, " Core state: ti=%3d bi=%3d energy=%8.5f "
705 "k=(%7.4f, %7.4f, %7.4f) rank=%d\n",
706 ti, bi, e, k[0], k[1], k[2], myComm->rank());
707 }
708 myComm->bcast (atom);
709 myComm->bcast(rmax);
710 myComm->bcast (l);
711 myComm->bcast (k);
712 int ng = g.size();
713 myComm->bcast(ng);
714 if (g.size() != ng)
715 {
716 g.resize(ng);
717 r.resize(ng);
718 }
719 myComm->bcast (g);
720 myComm->bcast (r);
721 double Z = (double)IonTypes[atom];
722 orbitalSet->MuffinTins[atom].addCore (l, m, r, g, k, Z);
723 icore++;
724 }
725 else
726 {
727 if (root)
728 {
729 int ti = SortBands[iorb].TwistIndex;
730 int bi = SortBands[iorb].BandIndex;
731 double e = SortBands[iorb].Energy;
732 PosType twist, k;
733 twist = TwistAngles[ti];
734 k = orbitalSet->PrimLattice.k_cart(twist);
735 fprintf (stderr, " Valence state: ti=%3d bi=%3d energy=%8.5f k=(%7.4f, %7.4f, %7.4f) rank=%d\n",
736 ti, bi, e, k[0], k[1], k[2], myComm->rank());
737 std::string vectorName = OrbitalPath (ti, bi) + "eigenvector";
738 HDFAttribIO<Array<std::complex<double>,3> > h_rawData(rawData);
739 h_rawData.read(H5FileID, vectorName.c_str());
740 if ((rawData.size(0) != nx) ||
741 (rawData.size(1) != ny) ||
742 (rawData.size(2) != nz))
743 {
744 fprintf (stderr, "Error in EinsplineSetBuilder::ReadBands.\n");
745 fprintf (stderr, "Extended orbitals should all have the same dimensions\n");
746 abort();
747 }
748 PosType ru;
749 for (int ix=0; ix<(nx-1); ix++)
750 {
751 ru[0] = (RealType)ix / (RealType)(nx-1);
752 for (int iy=0; iy<(ny-1); iy++)
753 {
754 ru[1] = (RealType)iy / (RealType)(ny-1);
755 for (int iz=0; iz<(nz-1); iz++)
756 {
757 ru[2] = (RealType)iz / (RealType)(nz-1);
758 double phi = -2.0*M_PI*dot (ru, TwistAngles[ti]);
759 double s, c;
760 qmcplusplus::sincos(phi, &s, &c);
761 std::complex<double> phase(c,s);
762 std::complex<double> z = phase*rawData(ix,iy,iz);
763 splineData(ix,iy,iz) = z.real();
764 }
765 }
766 }
767 }
768 myComm->bcast(splineData);
769 set_multi_UBspline_3d_d
770 (orbitalSet->MultiSpline, ival, splineData.data());
771 // Now read muffin tin data
772 for (int tin=0; tin<NumMuffinTins; tin++)
773 {
774 // app_log() << "Reading data for muffin tin " << tin << std::endl;
775 PosType twist, k;
776 int lmax = MT_APW_lmax[tin];
777 int numYlm = (lmax+1)*(lmax+1);
778 Array<std::complex<double>,2>
779 u_lm_r(numYlm, MT_APW_num_radial_points[tin]);
780 Array<std::complex<double>,1> du_lm_dr (numYlm);
781 if (root)
782 {
783 int ti = SortBands[iorb].TwistIndex;
784 int bi = SortBands[iorb].BandIndex;
785 twist = TwistAngles[ti];
786 k = orbitalSet->PrimLattice.k_cart(twist);
787 std::string uName = MuffinTinPath (ti, bi,tin) + "u_lm_r";
788 std::string duName = MuffinTinPath (ti, bi,tin) + "du_lm_dr";
789 HDFAttribIO<Array<std::complex<double>,2> > h_u_lm_r(u_lm_r);
790 HDFAttribIO<Array<std::complex<double>,1> > h_du_lm_dr(du_lm_dr);
791 h_u_lm_r.read(H5FileID, uName.c_str());
792 h_du_lm_dr.read(H5FileID, duName.c_str());
793 }
794 myComm->bcast(u_lm_r);
795 myComm->bcast(du_lm_dr);
796 myComm->bcast(k);
797 double Z = (double)IonTypes[tin];
798 orbitalSet->MuffinTins[tin].set_APW (ival, k, u_lm_r, du_lm_dr, Z);
799 }
800 ival++;
801 } // valence state
802 iorb++;
803 }
804 //ExtendedMap_d[set] = orbitalSet->MultiSpline;
805 }
806 #endif
807
808 } // namespace qmcplusplus
809