From 3dce32e450b6c181a0526f7fc7999a4ed893921f Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 21:31:10 -0500 Subject: [PATCH 1/9] Review EinsplineSetBuilder::AnalyzeTwists2 --- .../EinsplineSetBuilderCommon.cpp | 57 ++++++++++--------- 1 file changed, 31 insertions(+), 26 deletions(-) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index 2572243a6e..43ef9182e7 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -86,10 +86,7 @@ inline TinyVector FracPart(const TinyVector& twist) } -EinsplineSetBuilder::~EinsplineSetBuilder() -{ - DEBUG_MEMORY("EinsplineSetBuilder::~EinsplineSetBuilder"); -} +EinsplineSetBuilder::~EinsplineSetBuilder() { DEBUG_MEMORY("EinsplineSetBuilder::~EinsplineSetBuilder"); } bool EinsplineSetBuilder::CheckLattice() @@ -405,25 +402,30 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) S(i, j) = (double)TileMatrix(i, j); - std::vector superFracs; - std::vector> superSets; - { // build super twists - // This holds to which supercell kpoint each primitive k-point belongs - std::vector superIndex; - const int numPrimTwists = TwistAngles.size(); - for (int ki = 0; ki < numPrimTwists; ki++) + + const int num_prim_kpoints = TwistAngles.size(); + + // build a list of unique super twists that all the primitive cell k-point correspond to. + std::vector superFracs; // twist super twist coordinates + std::vector + superIndex; // the indices of the super twists that correpsond to all the primitive cell k-points in the unique list. + { + // scan all the primitive cell k-points + for (int ki = 0; ki < num_prim_kpoints; ki++) { PosType primTwist = TwistAngles[ki]; PosType superTwist = dot(S, primTwist); PosType kp = PrimCell.k_cart(primTwist); PosType ks = SuperCell.k_cart(superTwist); + // check the consistency of tiling, primitive and super cells. if (dot(ks - kp, ks - kp) > 1.0e-6) { app_error() << "Primitive and super k-points do not agree. Error in coding.\n"; APP_ABORT("EinsplineSetBuilder::AnalyzeTwists2"); } PosType frac = FracPart(superTwist); - bool found = false; + // verify if the super twist that correpsonds to this primitive cell k-point exists in the unique list or not. + bool found = false; for (int j = 0; j < superFracs.size(); j++) { PosType diff = frac - superFracs[j]; @@ -439,18 +441,14 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect superFracs.push_back(frac); } } - const int numSuperTwists = superFracs.size(); - app_log() << "Found " << numSuperTwists << " distinct supercell twists.\n"; - // For each supercell twist, create a list of primitive twists which - // belong to it. - superSets.resize(numSuperTwists); - for (int ki = 0; ki < numPrimTwists; ki++) - superSets[superIndex[ki]].push_back(ki); - app_log() << "number of things" << std::endl; - app_log() << TwistSymmetry.size() << std::endl; - app_log() << TwistWeight.size() << std::endl; - // for (int ki=0; ki 1 ? "s" : "") + << " based on " << num_prim_kpoints << " primitive cell k-point" << (num_prim_kpoints > 1 ? "s" : "") + << std::endl; if (myComm->rank() == 0) { int n_tot_irred(0); @@ -466,9 +464,16 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect } } } - const int numSuperTwists = superFracs.size(); - { // determine twist_num_ + // For each supercell twist, create a list of primitive twists which correspond to it. + std::vector> superSets; + { + superSets.resize(numSuperTwists); + for (int ki = 0; ki < num_prim_kpoints; ki++) + superSets[superIndex[ki]].push_back(ki); + } + + { // look up a super cell twist and return its index in the unique list of super cell twists. std::function find_twist = [&](const TinyVector& twist) { int twist_num = -1; PosType gtFrac = FracPart(twist); From 6b7570ddadf752d78edfe4d6857d8a0241944e1f Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 21:50:03 -0500 Subject: [PATCH 2/9] TwistAngle -> primcell_kpoints. --- .../BsplineFactory/BsplineReaderBase.cpp | 6 ++-- .../BsplineFactory/BsplineReaderBase.h | 4 +-- .../BsplineFactory/SplineSetReader.h | 4 +-- src/QMCWaveFunctions/EinsplineSetBuilder.h | 3 +- .../EinsplineSetBuilderCommon.cpp | 28 +++++++++---------- .../EinsplineSetBuilderESHDF.fft.cpp | 8 +++--- .../EinsplineSetBuilderOld.cpp | 10 +++---- .../EinsplineSetBuilder_createSPOs.cpp | 2 +- 8 files changed, 33 insertions(+), 32 deletions(-) diff --git a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.cpp b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.cpp index 2395d275c2..a387eab57f 100644 --- a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.cpp +++ b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.cpp @@ -220,10 +220,10 @@ void BsplineReaderBase::initialize_spo2band(int spin, int bi = bigspace[i].BandIndex; double e = bigspace[i].Energy; int nd = (bigspace[i].MakeTwoCopies) ? 2 : 1; - PosType k = mybuilder->PrimCell.k_cart(mybuilder->TwistAngles[ti]); + PosType k = mybuilder->PrimCell.k_cart(mybuilder->primcell_kpoints[ti]); int s_size = std::snprintf(s.data(), s.size(), "%8d %8d %8d %8d %12.6f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", - i, ns, ti, bi, e, k[0], k[1], k[2], mybuilder->TwistAngles[ti][0], - mybuilder->TwistAngles[ti][1], mybuilder->TwistAngles[ti][2], nd); + i, ns, ti, bi, e, k[0], k[1], k[2], mybuilder->primcell_kpoints[ti][0], + mybuilder->primcell_kpoints[ti][1], mybuilder->primcell_kpoints[ti][2], nd); if (s_size < 0) throw std::runtime_error("Error generating bandinfo"); o << s.data(); diff --git a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h index a41fd60445..53fe185baa 100644 --- a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h +++ b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h @@ -113,7 +113,7 @@ struct BsplineReaderBase for (int iorb = 0; iorb < N; iorb++) { int ti = cur_bands[iorb].TwistIndex; - bspline->kPoints[iorb] = mybuilder->PrimCell.k_cart(-mybuilder->TwistAngles[ti]); + bspline->kPoints[iorb] = mybuilder->PrimCell.k_cart(-mybuilder->primcell_kpoints[ti]); bspline->MakeTwoCopies[iorb] = (num < (numOrbs - 1)) && cur_bands[iorb].MakeTwoCopies; num += bspline->MakeTwoCopies[iorb] ? 2 : 1; } @@ -125,7 +125,7 @@ struct BsplineReaderBase if (!bspline->isComplex()) { //no k-point folding, single special k point (G, L ...) - TinyVector twist0 = mybuilder->TwistAngles[bandgroup.TwistIndex]; + TinyVector twist0 = mybuilder->primcell_kpoints[bandgroup.TwistIndex]; for (int i = 0; i < 3; i++) if (bconds[i] && ((std::abs(std::abs(twist0[i]) - 0.5) < 1.0e-8))) bspline->HalfG[i] = 1; diff --git a/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h b/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h index 7aedcf49b3..5d6153f7fa 100644 --- a/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h +++ b/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h @@ -259,7 +259,7 @@ struct SplineSetReader : public BsplineReaderBase if (bspline->isComplex()) { if (rotate) - fix_phase_rotate_c2c(FFTbox, splineData_r, splineData_i, mybuilder->TwistAngles[ti], rotate_phase_r, + fix_phase_rotate_c2c(FFTbox, splineData_r, splineData_i, mybuilder->primcell_kpoints[ti], rotate_phase_r, rotate_phase_i); else { @@ -272,7 +272,7 @@ struct SplineSetReader : public BsplineReaderBase } else { - fix_phase_rotate_c2r(FFTbox, splineData_r, mybuilder->TwistAngles[ti], rotate_phase_r, rotate_phase_i); + fix_phase_rotate_c2r(FFTbox, splineData_r, mybuilder->primcell_kpoints[ti], rotate_phase_r, rotate_phase_i); einspline::set(spline_r, splineData_r.data()); } } diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.h b/src/QMCWaveFunctions/EinsplineSetBuilder.h index f7e66293d3..ffcd9ccb7c 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder.h +++ b/src/QMCWaveFunctions/EinsplineSetBuilder.h @@ -214,7 +214,8 @@ class EinsplineSetBuilder : public SPOSetBuilder ///////////////////////////// // The "true" twist number after analyzing twistnum, twist XML input and h5 int twist_num_; - std::vector> TwistAngles; + // primitive cell k-points from DFT calculations + std::vector> primcell_kpoints; // integer index of sym operation from the irreducible brillion zone std::vector TwistSymmetry; // number of twists equivalent to this one in the big DFT grid diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index 43ef9182e7..30d20a32d1 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -196,10 +196,10 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() //myComm->bcast(IonTypes); bbuffer.add(&IonPos[0][0], &IonPos[0][0] + OHMMS_DIM * numIons); //myComm->bcast(IonPos); - if (TwistAngles.size() != NumTwists) - TwistAngles.resize(NumTwists); - bbuffer.add(&TwistAngles[0][0], &TwistAngles[0][0] + OHMMS_DIM * NumTwists); - //myComm->bcast(TwistAngles); + if (primcell_kpoints.size() != NumTwists) + primcell_kpoints.resize(NumTwists); + bbuffer.add(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); + //myComm->bcast(primcell_kpoints); if (TwistSymmetry.size() != NumTwists) TwistSymmetry.resize(NumTwists); bibuffer.add(&TwistSymmetry[0], &TwistSymmetry[0] + NumTwists); @@ -236,7 +236,7 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() for (int i = 0; i < numIons; ++i) bibuffer.get(IonTypes[i]); bbuffer.get(&IonPos[0][0], &IonPos[0][0] + OHMMS_DIM * numIons); - bbuffer.get(&TwistAngles[0][0], &TwistAngles[0][0] + OHMMS_DIM * NumTwists); + bbuffer.get(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); bibuffer.get(&TwistSymmetry[0], &TwistSymmetry[0] + NumTwists); bibuffer.get(&TwistWeight[0], &TwistWeight[0] + NumTwists); bbuffer.get(MT_APW_radii.begin(), MT_APW_radii.end()); @@ -403,7 +403,7 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect for (int j = 0; j < 3; j++) S(i, j) = (double)TileMatrix(i, j); - const int num_prim_kpoints = TwistAngles.size(); + const int num_prim_kpoints = primcell_kpoints.size(); // build a list of unique super twists that all the primitive cell k-point correspond to. std::vector superFracs; // twist super twist coordinates @@ -413,7 +413,7 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect // scan all the primitive cell k-points for (int ki = 0; ki < num_prim_kpoints; ki++) { - PosType primTwist = TwistAngles[ki]; + PosType primTwist = primcell_kpoints[ki]; PosType superTwist = dot(S, primTwist); PosType kp = PrimCell.k_cart(primTwist); PosType ks = SuperCell.k_cart(superTwist); @@ -577,12 +577,12 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect int N = superSets[si].size(); for (int i = 0; i < N; i++) { - PosType twistPrim_i = TwistAngles[superSets[si][i]]; + PosType twistPrim_i = primcell_kpoints[superSets[si][i]]; PosType twistSuper_i = dot(S, twistPrim_i); PosType superInt_i = IntPart(twistSuper_i); for (int j = i + 1; j < N; j++) { - PosType twistPrim_j = TwistAngles[superSets[si][j]]; + PosType twistPrim_j = primcell_kpoints[superSets[si][j]]; PosType twistSuper_j = dot(S, twistPrim_j); PosType superInt_j = IntPart(twistSuper_j); if (dot(superInt_i - superInt_j, superInt_i - superInt_j) < 1.0e-6) @@ -606,12 +606,12 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect for (int i = 0; i < IncludeTwists.size(); i++) { int ti = IncludeTwists[i]; - PosType twist_i = TwistAngles[ti]; + PosType twist_i = primcell_kpoints[ti]; bool distinct = true; for (int j = i + 1; j < IncludeTwists.size(); j++) { int tj = IncludeTwists[j]; - PosType twist_j = TwistAngles[tj]; + PosType twist_j = primcell_kpoints[tj]; PosType sum = twist_i + twist_j; PosType diff = twist_i - twist_j; if (TwistPair(twist_i, twist_j)) @@ -628,11 +628,11 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect { MakeTwoCopies[i] = false; int ti = DistinctTwists[i]; - PosType twist_i = TwistAngles[ti]; + PosType twist_i = primcell_kpoints[ti]; for (int j = 0; j < copyTwists.size(); j++) { int tj = copyTwists[j]; - PosType twist_j = TwistAngles[tj]; + PosType twist_j = primcell_kpoints[tj]; if (TwistPair(twist_i, twist_j)) MakeTwoCopies[i] = true; } @@ -652,7 +652,7 @@ void EinsplineSetBuilder::AnalyzeTwists2(const int twist_num_inp, const TinyVect for (int i = 0; i < DistinctTwists.size(); i++) { int ti = DistinctTwists[i]; - PosType twist = TwistAngles[ti]; + PosType twist = primcell_kpoints[ti]; for (int j = 0; j < OHMMS_DIM; j++) if (std::abs(twist[j] - 0.0) > MatchingTol && std::abs(twist[j] - 0.5) > MatchingTol && std::abs(twist[j] + 0.5) > MatchingTol) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp index cccae49b56..32c61a674d 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp @@ -264,16 +264,16 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) /////////////////////////// // Read the twist angles // /////////////////////////// - TwistAngles.resize(NumTwists); + primcell_kpoints.resize(NumTwists); TwistSymmetry.resize(NumTwists); TwistWeight.resize(NumTwists); for (int ti = 0; ti < NumTwists; ti++) { std::ostringstream path; path << "/electrons/kpoint_" << ti << "/reduced_k"; - TinyVector TwistAngles_DP; - H5File.read(TwistAngles_DP, path.str()); - TwistAngles[ti] = TwistAngles_DP; + TinyVector primcell_kpoints_DP; + H5File.read(primcell_kpoints_DP, path.str()); + primcell_kpoints[ti] = primcell_kpoints_DP; if ((Version[0] >= 2) and (Version[1] >= 1)) { std::ostringstream sym_path; diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp index 862dfaccc8..c4d316445a 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp @@ -138,7 +138,7 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) /////////////////////////// // Read the twist angles // /////////////////////////// - TwistAngles.resize(NumTwists); + primcell_kpoints.resize(NumTwists); for (int ti = 0; ti < NumTwists; ti++) { std::ostringstream path; @@ -146,11 +146,11 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) path << eigenstatesGroup << "/twist_" << ti << "/twist_angle"; else path << eigenstatesGroup << "/twist/twist_angle"; - TinyVector TwistAngles_DP; - H5File.read(TwistAngles_DP, path.str()); - TwistAngles[ti] = TwistAngles_DP; + TinyVector primcell_kpoints_DP; + H5File.read(primcell_kpoints_DP, path.str()); + primcell_kpoints[ti] = primcell_kpoints_DP; int length = std::snprintf(buff.data(), buff.size(), " Found twist angle (%6.3f, %6.3f, %6.3f)\n", - TwistAngles[ti][0], TwistAngles[ti][1], TwistAngles[ti][2]); + primcell_kpoints[ti][0], primcell_kpoints[ti][1], primcell_kpoints[ti][2]); if (length < 0) throw std::runtime_error("Error converting twist angle to string"); app_log() << std::string_view(buff.data(), length); diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp index 7e670592e6..223e689f58 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp @@ -346,7 +346,7 @@ std::unique_ptr EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur) for (int iorb = 0, num = 0; iorb < NumDistinctOrbitals; iorb++) { int ti = (*FullBands[spinSet])[iorb].TwistIndex; - temp_OrbitalSet->kPoints[iorb] = PrimCell.k_cart(-TwistAngles[ti]); + temp_OrbitalSet->kPoints[iorb] = PrimCell.k_cart(-primcell_kpoints[ti]); temp_OrbitalSet->MakeTwoCopies[iorb] = (num < (numOrbs - 1)) && (*FullBands[spinSet])[iorb].MakeTwoCopies; num += temp_OrbitalSet->MakeTwoCopies[iorb] ? 2 : 1; } From ed56bdc4db62805f7897c5f0e78632ebeeff5a3e Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 21:57:18 -0500 Subject: [PATCH 3/9] Remove TwistSymmetry TwistWeight --- src/QMCWaveFunctions/EinsplineSetBuilder.h | 4 ---- src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp | 9 --------- src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp | 11 ----------- 3 files changed, 24 deletions(-) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.h b/src/QMCWaveFunctions/EinsplineSetBuilder.h index ffcd9ccb7c..edfaf9fc89 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder.h +++ b/src/QMCWaveFunctions/EinsplineSetBuilder.h @@ -216,10 +216,6 @@ class EinsplineSetBuilder : public SPOSetBuilder int twist_num_; // primitive cell k-points from DFT calculations std::vector> primcell_kpoints; - // integer index of sym operation from the irreducible brillion zone - std::vector TwistSymmetry; - // number of twists equivalent to this one in the big DFT grid - std::vector TwistWeight; TinyVector TileFactor; Tensor TileMatrix; diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index 30d20a32d1..f35ffe73aa 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -199,13 +199,6 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() if (primcell_kpoints.size() != NumTwists) primcell_kpoints.resize(NumTwists); bbuffer.add(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); - //myComm->bcast(primcell_kpoints); - if (TwistSymmetry.size() != NumTwists) - TwistSymmetry.resize(NumTwists); - bibuffer.add(&TwistSymmetry[0], &TwistSymmetry[0] + NumTwists); - if (TwistWeight.size() != NumTwists) - TwistWeight.resize(NumTwists); - bibuffer.add(&TwistWeight[0], &TwistWeight[0] + NumTwists); bbuffer.add(MT_APW_radii.begin(), MT_APW_radii.end()); bibuffer.add(MT_APW_lmax.begin(), MT_APW_lmax.end()); bibuffer.add(MT_APW_num_radial_points.begin(), MT_APW_num_radial_points.end()); @@ -237,8 +230,6 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() bibuffer.get(IonTypes[i]); bbuffer.get(&IonPos[0][0], &IonPos[0][0] + OHMMS_DIM * numIons); bbuffer.get(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); - bibuffer.get(&TwistSymmetry[0], &TwistSymmetry[0] + NumTwists); - bibuffer.get(&TwistWeight[0], &TwistWeight[0] + NumTwists); bbuffer.get(MT_APW_radii.begin(), MT_APW_radii.end()); bibuffer.get(MT_APW_lmax.begin(), MT_APW_lmax.end()); bibuffer.get(MT_APW_num_radial_points.begin(), MT_APW_num_radial_points.end()); diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp index 32c61a674d..32f547a50b 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp @@ -265,8 +265,6 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) // Read the twist angles // /////////////////////////// primcell_kpoints.resize(NumTwists); - TwistSymmetry.resize(NumTwists); - TwistWeight.resize(NumTwists); for (int ti = 0; ti < NumTwists; ti++) { std::ostringstream path; @@ -274,15 +272,6 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) TinyVector primcell_kpoints_DP; H5File.read(primcell_kpoints_DP, path.str()); primcell_kpoints[ti] = primcell_kpoints_DP; - if ((Version[0] >= 2) and (Version[1] >= 1)) - { - std::ostringstream sym_path; - sym_path << "/electrons/kpoint_" << ti << "/symgroup"; - H5File.readEntry(TwistSymmetry[ti], sym_path.str()); - std::ostringstream nsym_path; - nsym_path << "/electrons/kpoint_" << ti << "/numsym"; - H5File.readEntry(TwistWeight[ti], nsym_path.str()); - } } if (qmc_common.use_density) { From 302ba2d2dd927fa20a4f4549bbc01dd780837b9b Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 22:51:13 -0500 Subject: [PATCH 4/9] Remove EinsplineSet and AtomicOrbital, --- src/QMCWaveFunctions/AtomicOrbital.cpp | 152 -- src/QMCWaveFunctions/AtomicOrbital.h | 657 --------- src/QMCWaveFunctions/BandInfo.h | 4 +- .../BsplineFactory/BsplineReaderBase.h | 7 +- .../BsplineFactory/SplineSetReader.h | 57 - src/QMCWaveFunctions/CMakeLists.txt | 4 - src/QMCWaveFunctions/EinsplineSet.cpp | 1285 ----------------- src/QMCWaveFunctions/EinsplineSet.h | 359 ----- src/QMCWaveFunctions/EinsplineSetBuilder.h | 23 +- .../EinsplineSetBuilderCommon.cpp | 94 +- .../EinsplineSetBuilderESHDF.fft.cpp | 65 +- .../EinsplineSetBuilderOld.cpp | 72 +- .../EinsplineSetBuilder_createSPOs.cpp | 104 +- .../EinsplineSpinorSetBuilder.cpp | 4 +- 14 files changed, 17 insertions(+), 2870 deletions(-) delete mode 100644 src/QMCWaveFunctions/AtomicOrbital.cpp delete mode 100644 src/QMCWaveFunctions/AtomicOrbital.h delete mode 100644 src/QMCWaveFunctions/EinsplineSet.cpp delete mode 100644 src/QMCWaveFunctions/EinsplineSet.h diff --git a/src/QMCWaveFunctions/AtomicOrbital.cpp b/src/QMCWaveFunctions/AtomicOrbital.cpp deleted file mode 100644 index 2042461443..0000000000 --- a/src/QMCWaveFunctions/AtomicOrbital.cpp +++ /dev/null @@ -1,152 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. -// -// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory -// -// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -////////////////////////////////////////////////////////////////////////////////////// - - -#include "AtomicOrbital.h" - -namespace qmcplusplus -{ -template<> -void AtomicOrbital>::allocate() -{ - Numlm = (lMax + 1) * (lMax + 1); - YlmVec.resize(Numlm); - dYlm_dthetaVec.resize(Numlm); - dYlm_dphiVec.resize(Numlm); - ulmVec.resize(Numlm * NumBands); - dulmVec.resize(Numlm * NumBands); - d2ulmVec.resize(Numlm * NumBands); - PolyCoefs.resize(PolyOrder + 1, NumBands, Numlm); - BCtype_z bc; - bc.lCode = NATURAL; - bc.rCode = NATURAL; - Ugrid grid; - grid.start = 0.0; - grid.end = SplineRadius; - grid.num = SplinePoints; - // if (RadialSpline) destroy_Bspline (RadialSpline); - RadialSpline = create_multi_UBspline_1d_z(grid, bc, Numlm * NumBands); - TwistAngles.resize(NumBands); -} - -template<> -void AtomicOrbital::allocate() -{ - Numlm = (lMax + 1) * (lMax + 1); - YlmVec.resize(Numlm); - dYlm_dthetaVec.resize(Numlm); - dYlm_dphiVec.resize(Numlm); - ulmVec.resize(Numlm * NumBands); - dulmVec.resize(Numlm * NumBands); - d2ulmVec.resize(Numlm * NumBands); - PolyCoefs.resize(PolyOrder + 1, NumBands, Numlm); - BCtype_d bc; - bc.lCode = NATURAL; - bc.rCode = NATURAL; - Ugrid grid; - grid.start = 0.0; - grid.end = SplineRadius; - grid.num = SplinePoints; - RadialSpline = create_multi_UBspline_1d_d(grid, bc, Numlm * NumBands); - TwistAngles.resize(NumBands); -} - - -template<> -void AtomicOrbital>::set_band(int band, - Array, 2>& spline_data, - Array, 2>& poly_coefs, - PosType twist) -{ - std::vector> one_spline(SplinePoints); - for (int lm = 0; lm < Numlm; lm++) - { - int index = band * Numlm + lm; - for (int i = 0; i < SplinePoints; i++) - one_spline[i] = spline_data(i, lm); - set_multi_UBspline_1d_z(RadialSpline, index, &one_spline[0]); - for (int n = 0; n <= PolyOrder; n++) - PolyCoefs(n, band, lm) = poly_coefs(n, lm); - } - TwistAngles[band] = twist; -} - - -// Here, we convert the complex Ylm representation to the real Ylm representation -template<> -void AtomicOrbital::set_band(int band, - Array, 2>& spline_data, - Array, 2>& poly_coefs, - PosType twist) -{ - std::vector one_spline(SplinePoints); - for (int l = 0; l <= lMax; l++) - { - // Set spline for m=0 - for (int i = 0; i < SplinePoints; i++) - one_spline[i] = spline_data(i, l * (l + 1)).real(); - int index = band * Numlm + l * (l + 1); - set_multi_UBspline_1d_d(RadialSpline, index, &one_spline[0]); - // Set poly ofr m=0 - for (int n = 0; n <= PolyOrder; n++) - PolyCoefs(n, band, l * (l + 1)) = poly_coefs(n, l * (l + 1)).real(); - // Set spline and poly for |m| > 0 - double minus_1_to_m = -1.0; - for (int m = 1; m <= l; m++) - { - int lmp = l * (l + 1) + m; - int lmm = l * (l + 1) - m; - index = band * Numlm + lmp; - for (int i = 0; i < SplinePoints; i++) - one_spline[i] = (spline_data(i, lmp).real() + minus_1_to_m * spline_data(i, lmm).real()); - set_multi_UBspline_1d_d(RadialSpline, index, &one_spline[0]); - index = band * Numlm + lmm; - for (int i = 0; i < SplinePoints; i++) - one_spline[i] = (-spline_data(i, lmp).imag() + minus_1_to_m * spline_data(i, lmm).imag()); - set_multi_UBspline_1d_d(RadialSpline, index, &one_spline[0]); - for (int n = 0; n <= PolyOrder; n++) - { - PolyCoefs(n, band, lmp) = (poly_coefs(n, lmp).real() + minus_1_to_m * poly_coefs(n, lmm).real()); - PolyCoefs(n, band, lmm) = (-poly_coefs(n, lmp).imag() + minus_1_to_m * poly_coefs(n, lmm).imag()); - } - minus_1_to_m *= -1.0; - } - } - TwistAngles[band] = twist; - // AtomicOrbital > zorb; - // zorb.set_pos (Pos); - // zorb.set_lmax(lMax); - // zorb.set_cutoff(CutoffRadius); - // zorb.set_spline(SplineRadius, SplinePoints); - // zorb.set_polynomial (PolyRadius, PolyOrder); - // zorb.set_num_bands(NumBands); - // zorb.allocate(); - // zorb.set_band(band, spline_data, poly_coefs, twist); - // PosType dir(0.324, -0.8, 1.3); - // dir = (1.0/std::sqrt(dot(dir,dir)))*dir; - // std::ostringstream fname; - // fname << "TestAtomic_" << band << ".dat"; - // FILE *fout = fopen (fname.str().c_str(), "w"); - // Vector zval(NumBands), val(NumBands); - // Vector zlapl(NumBands), lapl(NumBands); - // Vector zgrad(NumBands), grad(NumBands); - // for (double u=-1.00001; u<=1.0; u+= 0.001) { - // PosType r = u*CutoffRadius * dir + Pos; - // zorb.evaluate(r, zval, zgrad, zlapl); - // evaluate(r, val, grad, lapl); - // fprintf (fout, "%12.8f %12.8f %12.8f %14.8e %14.8e\n", - // r[0], r[1], r[2], lapl[band], lapl[band]); - // } - // fclose (fout); -} -} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/AtomicOrbital.h b/src/QMCWaveFunctions/AtomicOrbital.h deleted file mode 100644 index 1465d1c98b..0000000000 --- a/src/QMCWaveFunctions/AtomicOrbital.h +++ /dev/null @@ -1,657 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. -// -// File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Ye Luo, yeluo@anl.gov, Argonne National Laboratory -// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory -// -// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -////////////////////////////////////////////////////////////////////////////////////// - - -#ifndef ATOMIC_ORBITAL_H -#define ATOMIC_ORBITAL_H - -#include "CPU/math.hpp" -#include "einspline/multi_bspline.h" -#include "QMCWaveFunctions/SPOSet.h" -#include "Lattice/CrystalLattice.h" -#include -#include "Utilities/TimerManager.h" - - -namespace qmcplusplus -{ -/****************************************************************** -// This is just a template trick to avoid template specialization // -// in AtomicOrbital. // -******************************************************************/ - -template -struct AtomicOrbitalTraits -{}; -template<> -struct AtomicOrbitalTraits -{ - using SplineType = multi_UBspline_1d_d; -}; -template<> -struct AtomicOrbitalTraits> -{ - using SplineType = multi_UBspline_1d_z; -}; - -inline void EinsplineMultiEval(multi_UBspline_1d_d* spline, double x, double* val) -{ - eval_multi_UBspline_1d_d(spline, x, val); -} -inline void EinsplineMultiEval(multi_UBspline_1d_z* spline, double x, std::complex* val) -{ - eval_multi_UBspline_1d_z(spline, x, val); -} -inline void EinsplineMultiEval(multi_UBspline_1d_d* spline, double x, double* val, double* grad, double* lapl) -{ - eval_multi_UBspline_1d_d_vgl(spline, x, val, grad, lapl); -} -inline void EinsplineMultiEval(multi_UBspline_1d_z* spline, - double x, - std::complex* val, - std::complex* grad, - std::complex* lapl) -{ - eval_multi_UBspline_1d_z_vgl(spline, x, val, grad, lapl); -} - - -template -class AtomicOrbital -{ -public: - using PosType = QMCTraits::PosType; - using RealType = QMCTraits::RealType; - using UnitCellType = CrystalLattice; - using RealValueVector = Vector; - using RealGradVector = Vector>; - using ComplexValueVector = Vector>; - using ComplexGradVector = Vector, OHMMS_DIM>>; - using RealHessVector = Vector>; - using ComplexHessVector = Vector, OHMMS_DIM>>; - using SplineType = typename AtomicOrbitalTraits::SplineType; - -private: - // Store in order - // Index = l*(l+1) + m. There are (lMax+1)^2 Ylm's - std::vector YlmVec, dYlm_dthetaVec, dYlm_dphiVec, ulmVec, dulmVec, d2ulmVec; - - SplineType* RadialSpline; - // The first index is n in r^n, the second is lm = l*(l+1)+m - Array PolyCoefs; - NewTimer &YlmTimer, &SplineTimer, &SumTimer; - RealType rmagLast; - std::vector TwistAngles; - -public: - PosType Pos; - RealType CutoffRadius, SplineRadius, PolyRadius; - int SplinePoints; - int PolyOrder; - int lMax, Numlm, NumBands; - UnitCellType Lattice; - - inline void set_pos(PosType pos) { Pos = pos; } - inline void set_lmax(int lmax) { lMax = lmax; } - inline void set_cutoff(RealType cutoff) { CutoffRadius = cutoff; } - inline void set_spline(RealType radius, int points) - { - SplineRadius = radius; - SplinePoints = points; - } - inline void set_polynomial(RealType radius, int order) - { - PolyRadius = radius; - PolyOrder = order; - } - inline void set_num_bands(int num_bands) { NumBands = num_bands; } - SplineType* get_radial_spline() { return RadialSpline; } - Array& get_poly_coefs() { return PolyCoefs; } - - inline void registerTimers() - { - YlmTimer.reset(); - SplineTimer.reset(); - } - - void allocate(); - - void set_band(int band, - Array, 2>& spline_data, - Array, 2>& poly_coefs, - PosType twist); - inline void CalcYlm(PosType rhat, - std::vector>& Ylm, - std::vector>& dYlm_dtheta, - std::vector>& dYlm_dphi); - - inline void CalcYlm(PosType rhat, - std::vector& Ylm, - std::vector& dYlm_dtheta, - std::vector& dYlm_dphi); - - inline bool evaluate(PosType r, ComplexValueVector& vals); - inline bool evaluate(PosType r, ComplexValueVector& val, ComplexGradVector& grad, ComplexValueVector& lapl); - inline bool evaluate(PosType r, ComplexValueVector& val, ComplexGradVector& grad, ComplexHessVector& lapl); - inline bool evaluate(PosType r, RealValueVector& vals); - inline bool evaluate(PosType r, RealValueVector& val, RealGradVector& grad, RealValueVector& lapl); - inline bool evaluate(PosType r, RealValueVector& val, RealGradVector& grad, RealHessVector& lapl); - - - AtomicOrbital() - : RadialSpline(NULL), - YlmTimer(createGlobalTimer("AtomicOrbital::CalcYlm")), - SplineTimer(createGlobalTimer("AtomicOrbital::1D spline")), - SumTimer(createGlobalTimer("AtomicOrbital::Summation")), - rmagLast(std::numeric_limits::max()) - { - // Nothing else for now - } -}; - - -template -inline bool AtomicOrbital::evaluate(PosType r, ComplexValueVector& vals) -{ - PosType dr = r - Pos; - PosType u = Lattice.toUnit(dr); - PosType img; - for (int i = 0; i < OHMMS_DIM; i++) - { - img[i] = round(u[i]); - u[i] -= img[i]; - } - dr = Lattice.toCart(u); - double r2 = dot(dr, dr); - if (r2 > CutoffRadius * CutoffRadius) - return false; - double rmag = std::sqrt(r2); - PosType rhat = (1.0 / rmag) * dr; - // Evaluate Ylm's - CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); - if (std::abs(rmag - rmagLast) > 1.0e-6) - { - // Evaluate radial functions - if (rmag > PolyRadius) - EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0])); - else - { - for (int index = 0; index < ulmVec.size(); index++) - ulmVec[index] = StorageType(); - double r2n = 1.0; - for (int n = 0; n <= PolyOrder; n++) - { - int index = 0; - for (int i = 0; i < vals.size(); i++) - for (int lm = 0; lm < Numlm; lm++) - ulmVec[index++] += r2n * PolyCoefs(n, i, lm); - r2n *= rmag; - } - } - rmagLast = rmag; - } - SumTimer.start(); - int index = 0; - for (int i = 0; i < vals.size(); i++) - { - vals[i] = std::complex(); - for (int lm = 0; lm < Numlm; lm++) - vals[i] += ulmVec[index++] * YlmVec[lm]; - double phase = -2.0 * M_PI * dot(TwistAngles[i], img); - // fprintf (stderr, "phase[%d] = %1.2f pi\n", i, phase/M_PI); - // fprintf (stderr, "img = [%f,%f,%f]\n", img[0], img[1], img[2]); - double s, c; - qmcplusplus::sincos(phase, &s, &c); - vals[i] *= std::complex(c, s); - } - SumTimer.stop(); - return true; -} - - -template -inline bool AtomicOrbital::evaluate(PosType r, RealValueVector& vals) -{ - PosType dr = r - Pos; - PosType u = Lattice.toUnit(dr); - PosType img; - for (int i = 0; i < OHMMS_DIM; i++) - { - img[i] = round(u[i]); - u[i] -= img[i]; - } - dr = Lattice.toCart(u); - double r2 = dot(dr, dr); - if (r2 > CutoffRadius * CutoffRadius) - return false; - double rmag = std::sqrt(r2); - PosType rhat = (1.0 / rmag) * dr; - // Evaluate Ylm's - CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); - if (std::abs(rmag - rmagLast) > 1.0e-6) - { - // Evaluate radial functions - if (rmag > PolyRadius) - { - SplineTimer.start(); - EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0])); - SplineTimer.stop(); - } - else - { - for (int index = 0; index < ulmVec.size(); index++) - ulmVec[index] = StorageType(); - double r2n = 1.0; - for (int n = 0; n <= PolyOrder; n++) - { - int index = 0; - for (int i = 0; i < vals.size(); i++) - for (int lm = 0; lm < Numlm; lm++) - ulmVec[index++] += r2n * PolyCoefs(n, i, lm); - r2n *= rmag; - } - } - rmagLast = rmag; - } - SumTimer.start(); - int index = 0; - for (int i = 0; i < vals.size(); i++) - { - vals[i] = 0.0; - StorageType tmp = 0.0; - for (int lm = 0; lm < Numlm; lm++, index++) - tmp += ulmVec[index] * YlmVec[lm]; - //vals[i] += real(ulmVec[index++] * YlmVec[lm]); - // vals[i] += (ulmVec[index].real() * YlmVec[lm].real() - - // ulmVec[index].imag() * YlmVec[lm].imag()); - double phase = -2.0 * M_PI * dot(TwistAngles[i], img); - double s, c; - qmcplusplus::sincos(phase, &s, &c); - vals[i] = real(std::complex(c, s) * tmp); - } - SumTimer.stop(); - return true; -} - -template -inline bool AtomicOrbital::evaluate(PosType r, - RealValueVector& vals, - RealGradVector& grads, - RealHessVector& hess) -{ - APP_ABORT(" AtomicOrbital::evaluate not implemented for Hess. \n"); - return true; -} - - -template -inline bool AtomicOrbital::evaluate(PosType r, - RealValueVector& vals, - RealGradVector& grads, - RealValueVector& lapl) -{ - PosType dr = r - Pos; - PosType u = Lattice.toUnit(dr); - PosType img; - for (int i = 0; i < OHMMS_DIM; i++) - { - img[i] = round(u[i]); - u[i] -= img[i]; - } - dr = Lattice.toCart(u); - double r2 = dot(dr, dr); - if (r2 > CutoffRadius * CutoffRadius) - return false; - double rmag = std::sqrt(r2); - double rInv = 1.0 / rmag; - PosType rhat = rInv * dr; - double costheta = rhat[2]; - double sintheta = std::sqrt(1.0 - costheta * costheta); - double cosphi = rhat[0] / sintheta; - double sinphi = rhat[1] / sintheta; - PosType thetahat = PosType(costheta * cosphi, costheta * sinphi, -sintheta); - PosType phihat = PosType(-sinphi, cosphi, 0.0); - // Evaluate Ylm's - CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); - // Evaluate radial functions - if (rmag > PolyRadius) - { - SplineTimer.start(); - EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0]), &(dulmVec[0]), &(d2ulmVec[0])); - SplineTimer.stop(); - } - else - { - for (int index = 0; index < ulmVec.size(); index++) - { - ulmVec[index] = StorageType(); - dulmVec[index] = StorageType(); - d2ulmVec[index] = StorageType(); - } - double r2n = 1.0, r2nm1 = 0.0, r2nm2 = 0.0; - double dn = 0.0; - double dnm1 = -1.0; - for (int n = 0; n <= PolyOrder; n++) - { - int index = 0; - for (int i = 0; i < vals.size(); i++) - for (int lm = 0; lm < Numlm; lm++, index++) - { - StorageType c = PolyCoefs(n, i, lm); - ulmVec[index] += r2n * c; - dulmVec[index] += dn * r2nm1 * c; - d2ulmVec[index] += dn * dnm1 * r2nm2 * c; - } - dn += 1.0; - dnm1 += 1.0; - r2nm2 = r2nm1; - r2nm1 = r2n; - r2n *= rmag; - } - } - SumTimer.start(); - int index = 0; - for (int i = 0; i < vals.size(); i++) - { - vals[i] = 0.0; - for (int j = 0; j < OHMMS_DIM; j++) - grads[i][j] = 0.0; - lapl[i] = 0.0; - // Compute e^{-i k.L} phase factor - double phase = -2.0 * M_PI * dot(TwistAngles[i], img); - double s, c; - qmcplusplus::sincos(phase, &s, &c); - std::complex e2mikr(c, s); - StorageType tmp_val, tmp_lapl, grad_rhat, grad_thetahat, grad_phihat; - tmp_val = tmp_lapl = grad_rhat = grad_thetahat = grad_phihat = StorageType(); - int lm = 0; - for (int l = 0; l <= lMax; l++) - for (int m = -l; m <= l; m++, lm++, index++) - { - std::complex im(0.0, (double)m); - tmp_val += ulmVec[index] * YlmVec[lm]; - grad_rhat += dulmVec[index] * YlmVec[lm]; - grad_thetahat += ulmVec[index] * rInv * dYlm_dthetaVec[lm]; - grad_phihat += (ulmVec[index] * dYlm_dphiVec[lm]) / (rmag * sintheta); - //grad_phihat += (ulmVec[index] * im *YlmVec[lm])/(rmag*sintheta); - tmp_lapl += YlmVec[lm] * - (-(double)(l * (l + 1)) * rInv * rInv * ulmVec[index] + d2ulmVec[index] + 2.0 * rInv * dulmVec[index]); - } - vals[i] = real(e2mikr * tmp_val); - lapl[i] = real(e2mikr * tmp_lapl); - grads[i] = (real(e2mikr * grad_rhat) * rhat + real(e2mikr * grad_thetahat) * thetahat + - real(e2mikr * grad_phihat) * phihat); - } - SumTimer.stop(); - rmagLast = rmag; - return true; -} - -template -inline bool AtomicOrbital::evaluate(PosType r, - ComplexValueVector& vals, - ComplexGradVector& grads, - ComplexHessVector& hess) -{ - APP_ABORT(" AtomicOrbital::evaluate not implemented for Hess. \n"); - return true; -} - -template -inline bool AtomicOrbital::evaluate(PosType r, - ComplexValueVector& vals, - ComplexGradVector& grads, - ComplexValueVector& lapl) -{ - PosType dr = r - Pos; - PosType u = Lattice.toUnit(dr); - PosType img; - for (int i = 0; i < OHMMS_DIM; i++) - { - img[i] = round(u[i]); - u[i] -= img[i]; - } - dr = Lattice.toCart(u); - double r2 = dot(dr, dr); - if (r2 > CutoffRadius * CutoffRadius) - return false; - double rmag = std::sqrt(r2); - double rInv = 1.0 / rmag; - PosType rhat = rInv * dr; - double costheta = rhat[2]; - double sintheta = std::sqrt(1.0 - costheta * costheta); - double cosphi = rhat[0] / sintheta; - double sinphi = rhat[1] / sintheta; - PosType thetahat = PosType(costheta * cosphi, costheta * sinphi, -sintheta); - PosType phihat = PosType(-sinphi, cosphi, 0.0); - // Evaluate Ylm's - CalcYlm(rhat, YlmVec, dYlm_dthetaVec, dYlm_dphiVec); - // Evaluate radial functions - if (rmag > PolyRadius) - { - SplineTimer.start(); - EinsplineMultiEval(RadialSpline, rmag, &(ulmVec[0]), &(dulmVec[0]), &(d2ulmVec[0])); - SplineTimer.stop(); - } - else - { - for (int index = 0; index < ulmVec.size(); index++) - { - ulmVec[index] = StorageType(); - dulmVec[index] = StorageType(); - d2ulmVec[index] = StorageType(); - } - double r2n = 1.0, r2nm1 = 0.0, r2nm2 = 0.0; - double dn = 0.0; - double dnm1 = -1.0; - for (int n = 0; n <= PolyOrder; n++) - { - int index = 0; - for (int i = 0; i < vals.size(); i++) - for (int lm = 0; lm < Numlm; lm++, index++) - { - StorageType c = PolyCoefs(n, i, lm); - ulmVec[index] += r2n * c; - dulmVec[index] += dn * r2nm1 * c; - d2ulmVec[index] += dn * dnm1 * r2nm2 * c; - } - dn += 1.0; - dnm1 += 1.0; - r2nm2 = r2nm1; - r2nm1 = r2n; - r2n *= rmag; - } - } - SumTimer.start(); - int index = 0; - for (int i = 0; i < vals.size(); i++) - { - vals[i] = 0.0; - for (int j = 0; j < OHMMS_DIM; j++) - grads[i][j] = 0.0; - lapl[i] = 0.0; - int lm = 0; - StorageType grad_rhat, grad_thetahat, grad_phihat; - // Compute e^{-i k.L} phase factor - double phase = -2.0 * M_PI * dot(TwistAngles[i], img); - double s, c; - qmcplusplus::sincos(phase, &s, &c); - std::complex e2mikr(c, s); - for (int l = 0; l <= lMax; l++) - for (int m = -l; m <= l; m++, lm++, index++) - { - std::complex im(0.0, (double)m); - vals[i] += ulmVec[index] * YlmVec[lm]; - grad_rhat += dulmVec[index] * YlmVec[lm]; - grad_thetahat += ulmVec[index] * rInv * dYlm_dthetaVec[lm]; - grad_phihat += (ulmVec[index] * im * YlmVec[lm]) / (rmag * sintheta); - lapl[i] += YlmVec[lm] * - (-(double)(l * (l + 1)) * rInv * rInv * ulmVec[index] + d2ulmVec[index] + 2.0 * rInv * dulmVec[index]); - } - vals[i] *= e2mikr; - lapl[i] *= e2mikr; - for (int j = 0; j < OHMMS_DIM; j++) - { - grads[i][j] = e2mikr * (grad_rhat * rhat[j] + grad_thetahat * thetahat[j] + grad_phihat * phihat[j]); - } - } - SumTimer.stop(); - rmagLast = rmag; - return true; -} - - -// Fast implementation -// See Geophys. J. Int. (1998) 135,pp.307-309 -template -inline void AtomicOrbital::CalcYlm(PosType rhat, - std::vector>& Ylm, - std::vector>& dYlm_dtheta, - std::vector>& dYlm_dphi) -{ - YlmTimer.start(); - const double fourPiInv = 0.0795774715459477; - double costheta = rhat[2]; - double sintheta = std::sqrt(1.0 - costheta * costheta); - double cottheta = costheta / sintheta; - double cosphi, sinphi; - cosphi = rhat[0] / sintheta; - sinphi = rhat[1] / sintheta; - std::complex e2iphi(cosphi, sinphi); - double lsign = 1.0; - double dl = 0.0; - std::vector XlmVec(2 * lMax + 1), dXlmVec(2 * lMax + 1); - for (int l = 0; l <= lMax; l++) - { - XlmVec[2 * l] = lsign; - dXlmVec[2 * l] = dl * cottheta * XlmVec[2 * l]; - XlmVec[0] = lsign * XlmVec[2 * l]; - dXlmVec[0] = lsign * dXlmVec[2 * l]; - double dm = dl; - double msign = lsign; - for (int m = l; m > 0; m--) - { - double tmp = std::sqrt((dl + dm) * (dl - dm + 1.0)); - XlmVec[l + m - 1] = -(dXlmVec[l + m] + dm * cottheta * XlmVec[l + m]) / tmp; - dXlmVec[l + m - 1] = (dm - 1.0) * cottheta * XlmVec[l + m - 1] + XlmVec[l + m] * tmp; - // Copy to negative m - XlmVec[l - (m - 1)] = -msign * XlmVec[l + m - 1]; - dXlmVec[l - (m - 1)] = -msign * dXlmVec[l + m - 1]; - msign *= -1.0; - dm -= 1.0; - } - double sum = 0.0; - for (int m = -l; m <= l; m++) - sum += XlmVec[l + m] * XlmVec[l + m]; - // Now, renormalize the Ylms for this l - double norm = std::sqrt((2.0 * dl + 1.0) * fourPiInv / sum); - for (int m = -l; m <= l; m++) - { - XlmVec[l + m] *= norm; - dXlmVec[l + m] *= norm; - } - // Multiply by azimuthal phase and store in Ylm - std::complex e2imphi(1.0, 0.0); - std::complex eye(0.0, 1.0); - for (int m = 0; m <= l; m++) - { - Ylm[l * (l + 1) + m] = XlmVec[l + m] * e2imphi; - Ylm[l * (l + 1) - m] = XlmVec[l - m] * qmcplusplus::conj(e2imphi); - dYlm_dphi[l * (l + 1) + m] = (double)m * eye * XlmVec[l + m] * e2imphi; - dYlm_dphi[l * (l + 1) - m] = -(double)m * eye * XlmVec[l - m] * qmcplusplus::conj(e2imphi); - dYlm_dtheta[l * (l + 1) + m] = dXlmVec[l + m] * e2imphi; - dYlm_dtheta[l * (l + 1) - m] = dXlmVec[l - m] * qmcplusplus::conj(e2imphi); - e2imphi *= e2iphi; - } - dl += 1.0; - lsign *= -1.0; - } - YlmTimer.stop(); -} - -// Fast implementation -// See Geophys. J. Int. (1998) 135,pp.307-309 -template -inline void AtomicOrbital::CalcYlm(PosType rhat, - std::vector& Ylm, - std::vector& dYlm_dtheta, - std::vector& dYlm_dphi) -{ - YlmTimer.start(); - const double fourPiInv = 0.0795774715459477; - double costheta = rhat[2]; - double sintheta = std::sqrt(1.0 - costheta * costheta); - double cottheta = costheta / sintheta; - double cosphi, sinphi; - cosphi = rhat[0] / sintheta; - sinphi = rhat[1] / sintheta; - std::complex e2iphi(cosphi, sinphi); - double lsign = 1.0; - double dl = 0.0; - std::vector XlmVec(2 * lMax + 1), dXlmVec(2 * lMax + 1); - for (int l = 0; l <= lMax; l++) - { - XlmVec[2 * l] = lsign; - dXlmVec[2 * l] = dl * cottheta * XlmVec[2 * l]; - XlmVec[0] = lsign * XlmVec[2 * l]; - dXlmVec[0] = lsign * dXlmVec[2 * l]; - double dm = dl; - double msign = lsign; - for (int m = l; m > 0; m--) - { - double tmp = std::sqrt((dl + dm) * (dl - dm + 1.0)); - XlmVec[l + m - 1] = -(dXlmVec[l + m] + dm * cottheta * XlmVec[l + m]) / tmp; - dXlmVec[l + m - 1] = (dm - 1.0) * cottheta * XlmVec[l + m - 1] + XlmVec[l + m] * tmp; - // Copy to negative m - XlmVec[l - (m - 1)] = -msign * XlmVec[l + m - 1]; - dXlmVec[l - (m - 1)] = -msign * dXlmVec[l + m - 1]; - msign *= -1.0; - dm -= 1.0; - } - double sum = 0.0; - for (int m = -l; m <= l; m++) - sum += XlmVec[l + m] * XlmVec[l + m]; - // Now, renormalize the Ylms for this l - double norm = std::sqrt((2.0 * dl + 1.0) * fourPiInv / sum); - for (int m = -l; m <= l; m++) - { - XlmVec[l + m] *= norm; - dXlmVec[l + m] *= norm; - } - // Multiply by azimuthal phase and store in Ylm - Ylm[l * (l + 1)] = XlmVec[l]; - dYlm_dphi[l * (l + 1)] = 0.0; - dYlm_dtheta[l * (l + 1)] = dXlmVec[l]; - std::complex e2imphi = e2iphi; - for (int m = 1; m <= l; m++) - { - Ylm[l * (l + 1) + m] = XlmVec[l + m] * e2imphi.real(); - Ylm[l * (l + 1) - m] = XlmVec[l + m] * e2imphi.imag(); - dYlm_dphi[l * (l + 1) + m] = -(double)m * XlmVec[l + m] * e2imphi.imag(); - dYlm_dphi[l * (l + 1) - m] = (double)m * XlmVec[l + m] * e2imphi.real(); - dYlm_dtheta[l * (l + 1) + m] = dXlmVec[l + m] * e2imphi.real(); - dYlm_dtheta[l * (l + 1) - m] = dXlmVec[l + m] * e2imphi.imag(); - e2imphi *= e2iphi; - } - dl += 1.0; - lsign *= -1.0; - } - YlmTimer.stop(); -} - - -} // namespace qmcplusplus -#endif diff --git a/src/QMCWaveFunctions/BandInfo.h b/src/QMCWaveFunctions/BandInfo.h index 0e8e8acd61..9cd070586a 100644 --- a/src/QMCWaveFunctions/BandInfo.h +++ b/src/QMCWaveFunctions/BandInfo.h @@ -37,11 +37,9 @@ struct BandInfo double Energy; /// This is true if we should make distinct copies represeninting a +k, -k pair bool MakeTwoCopies; - /// True if this state is a core state - bool IsCoreState; ///default constructor BandInfo() - : TwistIndex(0), BandIndex(-1), BandGroup(0), Spin(0), Energy(1e9), MakeTwoCopies(false), IsCoreState(false) + : TwistIndex(0), BandIndex(-1), BandGroup(0), Spin(0), Energy(1e9), MakeTwoCopies(false) {} /** operator to determine the order of any band diff --git a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h index 53fe185baa..2a69c08aaf 100644 --- a/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h +++ b/src/QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h @@ -19,8 +19,11 @@ */ #ifndef QMCPLUSPLUS_BSPLINE_READER_BASE_H #define QMCPLUSPLUS_BSPLINE_READER_BASE_H + #include "mpi/collectives.h" #include "mpi/point2point.h" +#include + namespace qmcplusplus { struct SPOSetInputInfo; @@ -189,10 +192,6 @@ struct BsplineReaderBase const std::vector& bigspace, SPOSetInfo& sposet, std::vector& band2spo); - - /** export the MultiSpline to the old class EinsplineSetExtended for the GPU calculation*/ - virtual std::unique_ptr export_MultiSplineComplexDouble() = 0; - virtual std::unique_ptr export_MultiSplineDouble() = 0; }; } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h b/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h index 5d6153f7fa..dbeb68ff3c 100644 --- a/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h +++ b/src/QMCWaveFunctions/BsplineFactory/SplineSetReader.h @@ -68,63 +68,6 @@ struct SplineSetReader : public BsplineReaderBase // transform cG to radial functions virtual void create_atomic_centers_Gspace(Vector>& cG, Communicate& band_group_comm, int iorb) {} - /** for exporting data from multi_UBspline_3d_d to multi_UBspline_3d_z - * This is only used by the legacy EinsplineSet class. To be deleted together with EinsplineSet. - */ - std::unique_ptr export_MultiSplineComplexDouble() override - { - Ugrid xyz_grid[3]; - BCtype_d xyz_bc_d[3]; - set_grid(bspline->HalfG, xyz_grid, xyz_bc_d); - - BCtype_z xyz_bc[3]; - for (int i = 0; i < 3; i++) - { - xyz_bc[i].lCode = xyz_bc_d[i].lCode; - xyz_bc[i].rCode = xyz_bc_d[i].rCode; - } - - const auto* source = (multi_UBspline_3d_d*)bspline->SplineInst->getSplinePtr(); - std::unique_ptr target; - target.reset(einspline::create(target.get(), xyz_grid, xyz_bc, source->num_splines / 2)); - - if (source->x_grid.num != target->x_grid.num || source->y_grid.num != target->y_grid.num || - source->z_grid.num != target->z_grid.num) - throw std::runtime_error("export_MultiSplineComplexDouble failed for inconsistent grid dimensions."); - - if (source->coefs_size != target->coefs_size * 2) - throw std::runtime_error("export_MultiSplineComplexDouble failed for inconsistent coefs_size."); - - std::copy_n(source->coefs, source->coefs_size, (double*)target->coefs); - - return target; - } - - /** for exporting data from multi_UBspline_3d_d to multi_UBspline_3d_z - * This is only used by the legacy EinsplineSet class. To be deleted together with EinsplineSet. - */ - std::unique_ptr export_MultiSplineDouble() override - { - Ugrid xyz_grid[3]; - BCtype_d xyz_bc[3]; - set_grid(bspline->HalfG, xyz_grid, xyz_bc); - - const auto* source = (multi_UBspline_3d_d*)bspline->SplineInst->getSplinePtr(); - std::unique_ptr target; - target.reset(einspline::create(target.get(), xyz_grid, xyz_bc, source->num_splines)); - - if (source->x_grid.num != target->x_grid.num || source->y_grid.num != target->y_grid.num || - source->z_grid.num != target->z_grid.num) - throw std::runtime_error("export_MultiSplineDouble failed for inconsistent grid dimensions."); - - if (source->coefs_size != target->coefs_size) - throw std::runtime_error("export_MultiSplineDouble failed for inconsistent coefs_size."); - - std::copy_n(source->coefs, source->coefs_size, target->coefs); - - return target; - } - std::unique_ptr create_spline_set(const std::string& my_name, int spin, const BandInfoGroup& bandgroup) override diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 8a0611a6f5..fe26571d8f 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -80,14 +80,10 @@ if(OHMMS_DIM MATCHES 3) endif(QMC_COMPLEX) if(HAVE_EINSPLINE) - if(NOT MIXED_PRECISION) - set(FERMION_SRCS ${FERMION_SRCS} EinsplineSet.cpp) - endif(NOT MIXED_PRECISION) set(FERMION_SRCS ${FERMION_SRCS} EinsplineSetBuilderCommon.cpp EinsplineSetBuilderOld.cpp - AtomicOrbital.cpp EinsplineSetBuilderReadBands_ESHDF.cpp EinsplineSetBuilderESHDF.fft.cpp EinsplineSetBuilder_createSPOs.cpp diff --git a/src/QMCWaveFunctions/EinsplineSet.cpp b/src/QMCWaveFunctions/EinsplineSet.cpp deleted file mode 100644 index 95bed41a33..0000000000 --- a/src/QMCWaveFunctions/EinsplineSet.cpp +++ /dev/null @@ -1,1285 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. -// -// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory -// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory -// -// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -////////////////////////////////////////////////////////////////////////////////////// - - -#include "CPU/e2iphi.h" -#include "EinsplineSet.h" -#include "einspline/multi_bspline.h" -#include "CPU/math.hpp" -#include "type_traits/ConvertToReal.h" - -namespace qmcplusplus -{ -template -inline void EinsplineSetExtended::computePhaseFactors(const TinyVector& r) -{ - APP_ABORT("EinsplineSetExtended::computePhaseFactors called"); - for (int i = 0; i < kPoints.size(); i++) - phase[i] = -dot(r, kPoints[i]); - eval_e2iphi(kPoints.size(), phase.data(), eikr.data()); - //eval_e2iphi(phase,eikr); - //#ifdef HAVE_MKL - // for (int i=0; i(c,s); - // } - //#endif -} - - -EinsplineSet::UnitCellType EinsplineSet::GetLattice() { return SuperLattice; } - -void EinsplineSet::resetSourceParticleSet(ParticleSet& ions) {} - -void EinsplineSet::setOrbitalSetSize(int norbs) { OrbitalSetSize = norbs; } - -// Real evaluation functions -inline void EinsplineMultiEval(multi_UBspline_3d_d* restrict spline, - const TinyVector& r, - Vector& psi) -{ - eval_multi_UBspline_3d_d(spline, r[0], r[1], r[2], psi.data()); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_d* restrict spline, TinyVector r, std::vector& psi) -{ - eval_multi_UBspline_3d_d(spline, r[0], r[1], r[2], &(psi[0])); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_d* restrict spline, - const TinyVector& r, - Vector& psi, - Vector>& grad) -{ - eval_multi_UBspline_3d_d_vg(spline, r[0], r[1], r[2], psi.data(), (double*)grad.data()); -} - - -inline void EinsplineMultiEval(multi_UBspline_3d_d* restrict spline, - const TinyVector& r, - Vector& psi, - Vector>& grad, - Vector>& hess) -{ - eval_multi_UBspline_3d_d_vgh(spline, r[0], r[1], r[2], psi.data(), (double*)grad.data(), (double*)hess.data()); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_d* restrict spline, - const TinyVector& r, - Vector& psi, - Vector>& grad, - Vector>& hess, - Vector, 3>>& gradhess) -{ - eval_multi_UBspline_3d_d_vghgh(spline, r[0], r[1], r[2], psi.data(), (double*)grad.data(), (double*)hess.data(), - (double*)gradhess.data()); -} - - -////////////////////////////////// -// Complex evaluation functions // -////////////////////////////////// -inline void EinsplineMultiEval(multi_UBspline_3d_z* restrict spline, - const TinyVector& r, - Vector>& psi) -{ - eval_multi_UBspline_3d_z(spline, r[0], r[1], r[2], psi.data()); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_z* restrict spline, - const TinyVector& r, - Vector>& psi, - Vector, 3>>& grad) -{ - eval_multi_UBspline_3d_z_vg(spline, r[0], r[1], r[2], psi.data(), (std::complex*)grad.data()); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_z* restrict spline, - const TinyVector& r, - Vector>& psi, - Vector, 3>>& grad, - Vector, 3>>& hess) -{ - eval_multi_UBspline_3d_z_vgh(spline, r[0], r[1], r[2], psi.data(), (std::complex*)grad.data(), - (std::complex*)hess.data()); -} - -inline void EinsplineMultiEval(multi_UBspline_3d_z* restrict spline, - const TinyVector& r, - Vector>& psi, - Vector, 3>>& grad, - Vector, 3>>& hess, - Vector, 3>, 3>>& gradhess) -{ - eval_multi_UBspline_3d_z_vghgh(spline, r[0], r[1], r[2], psi.data(), (std::complex*)grad.data(), - (std::complex*)hess.data(), (std::complex*)gradhess.data()); -} - -template -void EinsplineSetExtended::setOrbitalSetSize(int norbs) -{ - OrbitalSetSize = norbs; -} - -#if !defined(QMC_COMPLEX) -template -void EinsplineSetExtended::evaluateValue(const ParticleSet& P, int iat, RealValueVector& psi) -{ - ValueTimer.start(); - const PosType& r(P.activeR(iat)); - // Check atomic orbitals - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_); - if (inAtom) - break; - } - StorageValueVector& valVec = storage_value_vector_; - if (!inAtom) - { - PosType ru(PrimLattice.toUnit(r)); - for (int i = 0; i < OHMMS_DIM; i++) - ru[i] -= std::floor(ru[i]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, valVec); - EinsplineTimer.stop(); - // Add e^ikr phase to B-spline orbitals - for (int j = 0; j < NumValenceOrbs; j++) - { - PosType k = kPoints[j]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - valVec[j] *= e_mikr; - } - } - const int N = storage_value_vector_.size(); - int psiIndex = 0; - for (int j = 0; j < N; j++) - { - std::complex psi_val = storage_value_vector_[j]; - psi[psiIndex] = real(psi_val); - psiIndex++; - if (MakeTwoCopies[j]) - { - psi[psiIndex] = imag(psi_val); - psiIndex++; - } - } - ValueTimer.stop(); -} - - -// This is an explicit specialization of the above for real orbitals -// with a real return value, i.e. simulations at the gamma or L -// point. -template<> -void EinsplineSetExtended::evaluateValue(const ParticleSet& P, int iat, RealValueVector& psi) -{ - ValueTimer.start(); - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, psi); - if (inAtom) - break; - } - if (!inAtom) - { - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int i = 0; i < OHMMS_DIM; i++) - { - RealType img = std::floor(ru[i]); - ru[i] -= img; - sign += HalfG[i] * (int)img; - } - // Check atomic orbitals - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, psi); - EinsplineTimer.stop(); - if (sign & 1) - for (int j = 0; j < psi.size(); j++) - psi[j] *= -1.0; - } - ValueTimer.stop(); -} - - -// Value, gradient, and laplacian -template -void EinsplineSetExtended::evaluateVGL(const ParticleSet& P, - int iat, - RealValueVector& psi, - RealGradVector& dpsi, - RealValueVector& d2psi) -{ - VGLTimer.start(); - const PosType& r(P.activeR(iat)); - std::complex eye(0.0, 1.0); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_, storage_grad_vector_, storage_lapl_vector_); - if (inAtom) - break; - } - StorageValueVector& valVec = storage_value_vector_; - StorageGradVector& gradVec = storage_grad_vector_; - StorageValueVector& laplVec = storage_lapl_vector_; - // Finally, copy into output vectors - int psiIndex = 0; - const int N = storage_value_vector_.size(); - for (int j = 0; j < N; j++) - { - std::complex psi_val, psi_lapl; - TinyVector, OHMMS_DIM> psi_grad; - psi_val = storage_value_vector_[j]; - psi_grad = storage_grad_vector_[j]; - psi_lapl = storage_lapl_vector_[j]; - psi[psiIndex] = real(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi[psiIndex][n] = real(psi_grad[n]); - d2psi[psiIndex] = real(psi_lapl); - psiIndex++; - if (MakeTwoCopies[j]) - { - psi[psiIndex] = imag(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi[psiIndex][n] = imag(psi_grad[n]); - d2psi[psiIndex] = imag(psi_lapl); - psiIndex++; - } - } - VGLTimer.stop(); -} - -template<> -void EinsplineSetExtended::evaluateVGL(const ParticleSet& P, - int iat, - RealValueVector& psi, - RealGradVector& dpsi, - RealValueVector& d2psi) -{ - VGLTimer.start(); - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, psi, dpsi, d2psi); - if (inAtom) - break; - } - if (!inAtom) - { - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int i = 0; i < OHMMS_DIM; i++) - { - RealType img = std::floor(ru[i]); - ru[i] -= img; - sign += HalfG[i] * (int)img; - } - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, psi, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - if (sign & 1) - for (int j = 0; j < psi.size(); j++) - { - psi[j] *= -1.0; - storage_grad_vector_[j] *= -1.0; - storage_hess_vector_[j] *= -1.0; - } - for (int i = 0; i < psi.size(); i++) - { - dpsi[i] = dot(PrimLattice.G, storage_grad_vector_[i]); - d2psi[i] = trace(storage_hess_vector_[i], GGt); - } - } - VGLTimer.stop(); -} - - -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealValueMatrix& d2psi) -{ - std::complex eye(0.0, 1.0); - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_, storage_grad_vector_, storage_lapl_vector_); - if (inAtom) - break; - } - StorageValueVector& valVec = storage_value_vector_; - StorageGradVector& gradVec = storage_grad_vector_; - StorageValueVector& laplVec = storage_lapl_vector_; - // Finally, copy into output vectors - int psiIndex = 0; - const int N = storage_value_vector_.size(); - for (int j = 0; j < N; j++) - { - std::complex psi_val, psi_lapl; - TinyVector, OHMMS_DIM> psi_grad; - psi_val = storage_value_vector_[j]; - psi_grad = storage_grad_vector_[j]; - psi_lapl = storage_lapl_vector_[j]; - psi(i, psiIndex) = real(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = real(psi_grad[n]); - d2psi(i, psiIndex) = real(psi_lapl); - psiIndex++; - // if (psiIndex >= dpsi.cols()) { - // std::cerr << "Error: out of bounds writing in EinsplineSet::evalate.\n" - // << "psiIndex = " << psiIndex << " dpsi.cols() = " << dpsi.cols() << std::endl; - // } - if (MakeTwoCopies[j]) - { - psi(i, psiIndex) = imag(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = imag(psi_grad[n]); - d2psi(i, psiIndex) = imag(psi_lapl); - psiIndex++; - } - } - } - VGLMatTimer.stop(); -} - -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi) -{ - std::complex eye(0.0, 1.0); - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - if (inAtom) - break; - } - StorageValueVector& valVec = storage_value_vector_; - StorageGradVector& gradVec = storage_grad_vector_; - StorageHessVector& hessVec = storage_hess_vector_; - Tensor, OHMMS_DIM> tmphs; - // Finally, copy into output vectors - int psiIndex = 0; - const int N = storage_value_vector_.size(); - for (int j = 0; j < N; j++) - { - std::complex psi_val; - TinyVector, OHMMS_DIM> psi_grad; - psi_val = storage_value_vector_[j]; - psi_grad = storage_grad_vector_[j]; - tmphs = storage_hess_vector_[j]; - psi(i, psiIndex) = real(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = real(psi_grad[n]); - //d2psi(i,psiIndex) = real(psi_lapl); - // FIX FIX FIX - for (int n = 0; n < OHMMS_DIM * OHMMS_DIM; n++) - grad_grad_psi(i, psiIndex)[n] = real(tmphs(n)); - psiIndex++; - // if (psiIndex >= dpsi.cols()) { - // std::cerr << "Error: out of bounds writing in EinsplineSet::evalate.\n" - // << "psiIndex = " << psiIndex << " dpsi.cols() = " << dpsi.cols() << std::endl; - // } - if (MakeTwoCopies[j]) - { - psi(i, psiIndex) = imag(psi_val); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = imag(psi_grad[n]); - //d2psi(i,psiIndex) = imag(psi_lapl); - for (int n = 0; n < OHMMS_DIM * OHMMS_DIM; n++) - grad_grad_psi(i, psiIndex)[n] = imag(tmphs(n)); - psiIndex++; - } - } - } - VGLMatTimer.stop(); -} - -template -void EinsplineSetExtended::evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat, - RealGradMatrix& dpsi) -{ - if (hasIonDerivs()) - { - // Loop over dimensions - for (int dim = 0; dim < OHMMS_DIM; dim++) - { - // Loop over electrons - for (int iel = first, i = 0; iel < last; iel++, i++) - { - const PosType& r(P.activeR(iel)); - PosType ru(PrimLattice.toUnit(r)); - assert(FirstOrderSplines[iat][dim]); - EinsplineMultiEval(FirstOrderSplines[iat][dim], ru, storage_value_vector_); - int dpsiIndex = 0; - for (int j = 0; j < NumValenceOrbs; j++) - { - PosType k = kPoints[j]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - storage_value_vector_[j] *= e_mikr; - dpsi(i, dpsiIndex)[dim] = real(storage_value_vector_[j]); - dpsiIndex++; - if (MakeTwoCopies[j]) - { - dpsi(i, dpsiIndex)[dim] = imag(storage_value_vector_[j]); - dpsiIndex++; - } - } - } - } - for (int i = 0; i < (last - first); i++) - for (int j = 0; j < (last - first); j++) - dpsi(i, j) = dot(PrimLattice.G, dpsi(i, j)); - } -} - - -// Evaluate the gradient w.r.t. to ion iat of the gradient and -// laplacian of the orbitals w.r.t. the electrons -template -void EinsplineSetExtended::evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - RealGradMatrix& dphi, - RealHessMatrix& dgrad_phi, - RealGradMatrix& dlapl_phi) -{ - if (hasIonDerivs()) - { - std::complex eye(0.0, 1.0); - // Loop over dimensions - for (int dim = 0; dim < OHMMS_DIM; dim++) - { - // Loop over electrons - for (int iel = first, i = 0; iel < last; iel++, i++) - { - const PosType& r(P.activeR(iel)); - PosType ru(PrimLattice.toUnit(r)); - assert(FirstOrderSplines[iat_src][dim]); - EinsplineMultiEval(FirstOrderSplines[iat_src][dim], ru, storage_value_vector_, storage_grad_vector_, - storage_hess_vector_); - int dphiIndex = 0; - for (int j = 0; j < NumValenceOrbs; j++) - { - storage_grad_vector_[j] = dot(PrimLattice.G, storage_grad_vector_[j]); - storage_lapl_vector_[j] = trace(storage_hess_vector_[j], GGt); - std::complex u = storage_value_vector_[j]; - TinyVector, OHMMS_DIM> gradu = storage_grad_vector_[j]; - std::complex laplu = storage_lapl_vector_[j]; - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - storage_value_vector_[j] = e_mikr * u; - storage_grad_vector_[j] = e_mikr * (-eye * u * ck + gradu); - storage_lapl_vector_[j] = e_mikr * (-dot(k, k) * u - 2.0 * eye * dot(ck, gradu) + laplu); - dphi(i, dphiIndex)[dim] = real(storage_value_vector_[j]); - for (int k = 0; k < OHMMS_DIM; k++) - dgrad_phi(dphiIndex)[dim] = real(storage_grad_vector_[j][k]); - dlapl_phi(dphiIndex)[dim] = real(storage_lapl_vector_[j]); - dphiIndex++; - if (MakeTwoCopies[j]) - { - dphi(i, dphiIndex)[dim] = imag(storage_value_vector_[j]); - for (int k = 0; k < OHMMS_DIM; k++) - dgrad_phi(i, dphiIndex)(dim, k) = imag(storage_grad_vector_[j][k]); - dlapl_phi(i, dphiIndex)[dim] = imag(storage_lapl_vector_[j]); - dphiIndex++; - } - } - } - } - for (int i = 0; i < (last - first); i++) - for (int j = 0; j < (last - first); j++) - { - dphi(i, j) = dot(PrimLattice.G, dphi(i, j)); - // Check this one! - dgrad_phi(i, j) = dot(PrimLattice.G, dgrad_phi(i, j)); - dlapl_phi(i, j) = dot(PrimLattice.G, dlapl_phi(i, j)); - } - } -} - - -template<> -void EinsplineSetExtended::evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - RealGradMatrix& dphi, - RealHessMatrix& dgrad_phi, - RealGradMatrix& dlapl_phi) -{ - if (hasIonDerivs()) - { - // Loop over dimensions - for (int dim = 0; dim < OHMMS_DIM; dim++) - { - assert(FirstOrderSplines[iat_src][dim]); - // Loop over electrons - for (int iel = first, i = 0; iel < last; iel++, i++) - { - const PosType& r(P.activeR(iel)); - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int n = 0; n < OHMMS_DIM; n++) - { - RealType img = std::floor(ru[n]); - ru[n] -= img; - sign += HalfG[n] * (int)img; - } - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineMultiEval(FirstOrderSplines[iat_src][dim], ru, storage_value_vector_, storage_grad_vector_, - storage_hess_vector_); - if (sign & 1) - for (int j = 0; j < OrbitalSetSize; j++) - { - dphi(i, j)[dim] = -1.0 * storage_value_vector_[j]; - PosType g = -1.0 * dot(PrimLattice.G, storage_grad_vector_[j]); - for (int k = 0; k < OHMMS_DIM; k++) - dgrad_phi(i, j)(dim, k) = g[k]; - dlapl_phi(i, j)[dim] = -1.0 * trace(storage_hess_vector_[j], GGt); - } - else - for (int j = 0; j < OrbitalSetSize; j++) - { - dphi(i, j)[dim] = storage_value_vector_[j]; - PosType g = dot(PrimLattice.G, storage_grad_vector_[j]); - for (int k = 0; k < OHMMS_DIM; k++) - dgrad_phi(i, j)(dim, k) = g[k]; - dlapl_phi(i, j)[dim] = trace(storage_hess_vector_[j], GGt); - } - } - } - for (int i = 0; i < (last - first); i++) - for (int j = 0; j < (last - first); j++) - { - dphi(i, j) = dot(PrimLattice.G, dphi(i, j)); - // Check this one! - dgrad_phi(i, j) = dot(PrimLattice.G, dgrad_phi(i, j)); - dlapl_phi(i, j) = dot(PrimLattice.G, dlapl_phi(i, j)); - } - } -} -template<> -void EinsplineSetExtended::evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat, - RealGradMatrix& dpsi) -{ - if (hasIonDerivs()) - { - // Loop over dimensions - for (int dim = 0; dim < OHMMS_DIM; dim++) - { - assert(FirstOrderSplines[iat][dim]); - // Loop over electrons - for (int iel = first, i = 0; iel < last; iel++, i++) - { - const PosType& r(P.activeR(iel)); - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int n = 0; n < OHMMS_DIM; n++) - { - RealType img = std::floor(ru[n]); - ru[n] -= img; - sign += HalfG[n] * (int)img; - } - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineMultiEval(FirstOrderSplines[iat][dim], ru, storage_value_vector_); - if (sign & 1) - for (int j = 0; j < OrbitalSetSize; j++) - dpsi(i, j)[dim] = -1.0 * storage_value_vector_[j]; - else - for (int j = 0; j < OrbitalSetSize; j++) - dpsi(i, j)[dim] = storage_value_vector_[j]; - } - } - for (int i = 0; i < (last - first); i++) - for (int j = 0; j < (last - first); j++) - { - dpsi(i, j) = dot(PrimLattice.G, dpsi(i, j)); - } - } -} - -#else - -template -void EinsplineSetExtended::evaluateValue(const ParticleSet& P, int iat, ComplexValueVector& psi) -{ - ValueTimer.start(); - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int i = 0; i < OHMMS_DIM; i++) - ru[i] -= std::floor(ru[i]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_); - EinsplineTimer.stop(); - //computePhaseFactors(r); - for (int i = 0; i < psi.size(); i++) - { - PosType k = kPoints[i]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi[i] = e_mikr * storage_value_vector_[i]; - } - ValueTimer.stop(); -} - -// Value, gradient, and laplacian -template -void EinsplineSetExtended::evaluateVGL(const ParticleSet& P, - int iat, - ComplexValueVector& psi, - ComplexGradVector& dpsi, - ComplexValueVector& d2psi) -{ - VGLTimer.start(); - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int i = 0; i < OHMMS_DIM; i++) - ru[i] -= std::floor(ru[i]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - //computePhaseFactors(r); - std::complex eye(0.0, 1.0); - for (int j = 0; j < psi.size(); j++) - { - std::complex u, laplu; - TinyVector, OHMMS_DIM> gradu; - u = storage_value_vector_[j]; - gradu = dot(PrimLattice.G, storage_grad_vector_[j]); - laplu = trace(storage_hess_vector_[j], GGt); - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi[j] = e_mikr * u; - dpsi[j] = e_mikr * (-eye * u * ck + gradu); - //convertVec(e_mikr*(-eye*u*ck + gradu), dpsi[j]); - d2psi[j] = e_mikr * (-dot(k, k) * u - 2.0 * eye * dot(ck, gradu) + laplu); - } - VGLTimer.stop(); -} - -// Value, gradient, and laplacian -template -void EinsplineSetExtended::evaluateVGH(const ParticleSet& P, - int iat, - ComplexValueVector& psi, - ComplexGradVector& dpsi, - ComplexHessVector& grad_grad_psi) -{ - VGLTimer.start(); - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int i = 0; i < OHMMS_DIM; i++) - ru[i] -= std::floor(ru[i]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - //computePhaseFactors(r); - std::complex eye(0.0, 1.0); - for (int j = 0; j < psi.size(); j++) - { - std::complex u; - TinyVector, OHMMS_DIM> gradu; - Tensor, OHMMS_DIM> hs, tmphs; - u = storage_value_vector_[j]; - gradu = dot(PrimLattice.G, storage_grad_vector_[j]); - ////laplu = trace(storage_hess_vector_[j], GGt); - tmphs = dot(PrimLattice.G, storage_hess_vector_[j]); - //hs = dot(tmphs,PrimLattice.G); - hs = dot(tmphs, PrimLattice.Gt); - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi[j] = e_mikr * u; - dpsi[j] = e_mikr * (-eye * u * ck + gradu); - //convertVec(e_mikr*(-eye*u*ck + gradu), dpsi[j]); - //d2psi[j] = e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu); - grad_grad_psi[j] = - e_mikr * (hs - u * outerProduct(ck, ck) - eye * outerProduct(ck, gradu) - eye * outerProduct(gradu, ck)); - } - VGLTimer.stop(); -} -#endif - -#if !defined(QMC_COMPLEX) -template<> -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealValueMatrix& d2psi) -{ - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_, storage_grad_vector_, storage_lapl_vector_); - if (inAtom) - { - for (int j = 0; j < OrbitalSetSize; j++) - { - psi(i, j) = storage_value_vector_[j]; - dpsi(i, j) = storage_grad_vector_[j]; - d2psi(i, j) = storage_lapl_vector_[j]; - } - break; - } - } - if (!inAtom) - { - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int n = 0; n < OHMMS_DIM; n++) - { - RealType img = std::floor(ru[n]); - ru[n] -= img; - sign += HalfG[n] * (int)img; - } - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - if (sign & 1) - for (int j = 0; j < OrbitalSetSize; j++) - { - storage_value_vector_[j] *= -1.0; - storage_grad_vector_[j] *= -1.0; - storage_hess_vector_[j] *= -1.0; - } - for (int j = 0; j < OrbitalSetSize; j++) - { - psi(i, j) = storage_value_vector_[j]; - dpsi(i, j) = dot(PrimLattice.G, storage_grad_vector_[j]); - d2psi(i, j) = trace(storage_hess_vector_[j], GGt); - } - } - } - VGLMatTimer.stop(); -} - -template<> -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi) -{ - //APP_ABORT("evaluate_notranspose: Check Hessian, then remove this error message.\n") - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - bool inAtom = false; - for (int jat = 0; jat < AtomicOrbitals.size(); jat++) - { - inAtom = AtomicOrbitals[jat].evaluate(r, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - if (inAtom) - { - for (int j = 0; j < OrbitalSetSize; j++) - { - psi(i, j) = storage_value_vector_[j]; - dpsi(i, j) = storage_grad_vector_[j]; - grad_grad_psi(i, j) = storage_hess_vector_[j]; - } - break; - } - } - if (!inAtom) - { - PosType ru(PrimLattice.toUnit(r)); - int sign = 0; - for (int n = 0; n < OHMMS_DIM; n++) - { - RealType img = std::floor(ru[n]); - ru[n] -= img; - sign += HalfG[n] * (int)img; - } - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - if (sign & 1) - for (int j = 0; j < OrbitalSetSize; j++) - { - storage_value_vector_[j] *= -1.0; - storage_grad_vector_[j] *= -1.0; - storage_hess_vector_[j] *= -1.0; - } - for (int j = 0; j < OrbitalSetSize; j++) - { - psi(i, j) = storage_value_vector_[j]; - dpsi(i, j) = dot(PrimLattice.G, storage_grad_vector_[j]); - grad_grad_psi(i, j) = dot(PrimLattice.G, dot(storage_hess_vector_[j], PrimLattice.Gt)); - } - } - } - VGLMatTimer.stop(); -} - -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi, - RealGGGMatrix& grad_grad_grad_logdet) -{ - // APP_ABORT(" EinsplineSetExtended::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n"); - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_, - storage_grad_hess_vector_); - EinsplineTimer.stop(); - for (int j = 0; j < NumValenceOrbs; j++) - { - TinyVector, OHMMS_DIM> tmpg; - Tensor, OHMMS_DIM> tmphs; - TinyVector, OHMMS_DIM>, OHMMS_DIM> tmpghs; - tmpg = dot(PrimLattice.G, storage_grad_vector_[j]); - storage_grad_vector_[j] = tmpg; - tmphs = dot(PrimLattice.G, storage_hess_vector_[j]); - storage_hess_vector_[j] = dot(tmphs, PrimLattice.Gt); - for (int n = 0; n < OHMMS_DIM; n++) - { - tmphs = dot(PrimLattice.G, storage_grad_hess_vector_[j][n]); - tmpghs[n] = dot(tmphs, PrimLattice.Gt); - } - storage_grad_hess_vector_[j] = dot(PrimLattice.G, tmpghs); - } - std::complex eye(0.0, 1.0); - // StorageValueVector &valVec = - // storage_value_vector_; - // StorageGradVector &gradVec = - // storage_grad_vector_; - // StorageHessVector &hessVec = - // storage_hess_vector_; - // Tensor,OHMMS_DIM> tmphs; - for (int j = 0; j < NumValenceOrbs; j++) - { - // std::complex u = valVec[j]; - // TinyVector,OHMMS_DIM> gradu = gradVec[j]; - // tmphs = hessVec[j]; - // PosType k = kPoints[j]; - // TinyVector,OHMMS_DIM> ck; - // for (int n=0; n e_mikr (c,s); - // valVec[j] = e_mikr*u; - // gradVec[j] = e_mikr*(-eye*u*ck + gradu); - // hessVec[j] = e_mikr*(tmphs -u*outerProduct(ck,ck) - eye*outerProduct(ck,gradu) - eye*outerProduct(gradu,ck)); - std::complex u = (storage_value_vector_[j]); - TinyVector, OHMMS_DIM> gradu = (storage_grad_vector_[j]); - Tensor, OHMMS_DIM> tmphs = (storage_hess_vector_[j]); - // TinyVector,OHMMS_DIM>,OHMMS_DIM> tmpghs=(storage_grad_hess_vector_[j]); - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - storage_value_vector_[j] = e_mikr * u; - storage_grad_vector_[j] = e_mikr * (-eye * u * ck + gradu); - storage_hess_vector_[j] = - e_mikr * (tmphs - u * outerProduct(ck, ck) - eye * outerProduct(ck, gradu) - eye * outerProduct(gradu, ck)); - //Is this right? - storage_grad_hess_vector_[j] *= e_mikr; - for (unsigned a0(0); a0 < OHMMS_DIM; a0++) - for (unsigned a1(0); a1 < OHMMS_DIM; a1++) - for (unsigned a2(0); a2 < OHMMS_DIM; a2++) - storage_grad_hess_vector_[j][a0](a1, a2) += e_mikr * - (-1.0 * eye * (ck[a0] * tmphs(a1, a2) + ck[a1] * tmphs(a0, a2) + ck[a2] * tmphs(a0, a1)) - - (ck[a0] * ck[a1] * gradu[a2] + ck[a0] * ck[a2] * gradu[a1] + ck[a1] * ck[a2] * gradu[a0]) + - eye * ck[a0] * ck[a1] * ck[a2] * u); - } - int psiIndex(0); - for (int j = 0; j < NumValenceOrbs; j++) - { - if (MakeTwoCopies[j]) - { - psi(i, psiIndex) = imag(storage_value_vector_[j]); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = imag(storage_grad_vector_[j][n]); - for (int n = 0; n < OHMMS_DIM * OHMMS_DIM; n++) - grad_grad_psi(i, psiIndex)[n] = imag(storage_hess_vector_[j](n)); - for (int n = 0; n < OHMMS_DIM; n++) - for (int m = 0; m < OHMMS_DIM * OHMMS_DIM; m++) - grad_grad_grad_logdet(i, psiIndex)[n][m] = imag(storage_grad_hess_vector_[j][n](m)); - psiIndex++; - psi(i, psiIndex) = real(storage_value_vector_[j]); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = real(storage_grad_vector_[j][n]); - for (int n = 0; n < OHMMS_DIM * OHMMS_DIM; n++) - grad_grad_psi(i, psiIndex)[n] = real(storage_hess_vector_[j](n)); - for (int n = 0; n < OHMMS_DIM; n++) - for (int m = 0; m < OHMMS_DIM * OHMMS_DIM; m++) - grad_grad_grad_logdet(i, psiIndex)[n][m] = real(storage_grad_hess_vector_[j][n](m)); - psiIndex++; - } - else - { - psi(i, psiIndex) = real(storage_value_vector_[j]); - for (int n = 0; n < OHMMS_DIM; n++) - dpsi(i, psiIndex)[n] = real(storage_grad_vector_[j][n]); - for (int n = 0; n < OHMMS_DIM * OHMMS_DIM; n++) - grad_grad_psi(i, psiIndex)[n] = real(storage_hess_vector_[j](n)); - for (int n = 0; n < OHMMS_DIM; n++) - for (int m = 0; m < OHMMS_DIM * OHMMS_DIM; m++) - grad_grad_grad_logdet(i, psiIndex)[n][m] = real(storage_grad_hess_vector_[j][n](m)); - psiIndex++; - } - } - } - VGLMatTimer.stop(); -} - -template<> -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi, - RealGGGMatrix& grad_grad_grad_logdet) -{ - // APP_ABORT(" EinsplineSetExtended::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n"); - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - bool inAtom = false; - int psiIndex = 0; - const int N = storage_value_vector_.size(); - for (int j = 0; j < N; j++) - { - psi(i, psiIndex) = storage_value_vector_[j]; - dpsi(i, psiIndex) = dot(storage_grad_vector_[j], PrimLattice.G); - grad_grad_psi(i, psiIndex) = storage_hess_vector_[j]; - grad_grad_grad_logdet(i, psiIndex) = dot(storage_grad_hess_vector_[j], PrimLattice.G); - psiIndex++; - } - } - VGLMatTimer.stop(); -} - - -#else -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexValueMatrix& d2psi) -{ - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - //computePhaseFactors(r); - std::complex eye(0.0, 1.0); - for (int j = 0; j < psi.cols(); j++) - { - std::complex u, laplu; - TinyVector, OHMMS_DIM> gradu; - u = storage_value_vector_[j]; - gradu = dot(PrimLattice.G, storage_grad_vector_[j]); - laplu = trace(storage_hess_vector_[j], GGt); - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi(i, j) = e_mikr * u; - //psi(j,i) = e_mikr * u; - dpsi(i, j) = e_mikr * (-eye * u * ck + gradu); - //convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j)); - d2psi(i, j) = e_mikr * (-dot(k, k) * u - 2.0 * eye * dot(ck, gradu) + laplu); - } - } - VGLMatTimer.stop(); -} - -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexHessMatrix& grad_grad_psi) -{ - VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_); - EinsplineTimer.stop(); - //computePhaseFactors(r); - std::complex eye(0.0, 1.0); - for (int j = 0; j < OrbitalSetSize; j++) - { - std::complex u; - TinyVector, OHMMS_DIM> gradu; - Tensor, OHMMS_DIM> hs, tmphs; - u = storage_value_vector_[j]; - gradu = dot(PrimLattice.G, storage_grad_vector_[j]); - // tmphs = dot(transpose(PrimLattice.G),storage_hess_vector_[j]); - tmphs = dot(PrimLattice.G, storage_hess_vector_[j]); - hs = dot(tmphs, PrimLattice.Gt); - //laplu = trace(storage_hess_vector_[j], GGt); - PosType k = kPoints[j]; - TinyVector, OHMMS_DIM> ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi(i, j) = e_mikr * u; - //psi(j,i) = e_mikr * u; - dpsi(i, j) = e_mikr * (-eye * u * ck + gradu); - //convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j)); - grad_grad_psi(i, j) = - e_mikr * (hs - u * outerProduct(ck, ck) - eye * outerProduct(ck, gradu) - eye * outerProduct(gradu, ck)); - } - } - VGLMatTimer.stop(); -} - -template<> -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexHessMatrix& grad_grad_psi, - ComplexGGGMatrix& grad_grad_grad_logdet) -{ - APP_ABORT( - " EinsplineSetExtended::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n"); -} - -template -void EinsplineSetExtended::evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexHessMatrix& grad_grad_psi, - ComplexGGGMatrix& grad_grad_grad_logdet) -{ - // APP_ABORT(" EinsplineSetExtended::evaluate_notranspose not implemented for grad_grad_grad_logdet yet. \n"); - // VGLMatTimer.start(); - for (int iat = first, i = 0; iat < last; iat++, i++) - { - const PosType& r(P.activeR(iat)); - PosType ru(PrimLattice.toUnit(r)); - for (int n = 0; n < OHMMS_DIM; n++) - ru[n] -= std::floor(ru[n]); - EinsplineTimer.start(); - EinsplineMultiEval(MultiSpline, ru, storage_value_vector_, storage_grad_vector_, storage_hess_vector_, - storage_grad_hess_vector_); - EinsplineTimer.stop(); - Tensor PG; - PG = PrimLattice.G; - Tensor TPG; - TPG = transpose(PrimLattice.G); - Tensor, OHMMS_DIM> hs, tmphs; - TinyVector, OHMMS_DIM>, OHMMS_DIM> tmpghs, hvdot; - for (int j = 0; j < NumValenceOrbs; j++) - { - storage_grad_vector_[j] = dot(PG, storage_grad_vector_[j]); - tmphs = dot(PG, storage_hess_vector_[j]); - storage_hess_vector_[j] = dot(tmphs, TPG); - for (int n = 0; n < OHMMS_DIM; n++) - { - tmpghs[n] = dot(PG, storage_grad_hess_vector_[j][n]); - storage_grad_hess_vector_[j][n] = dot(tmpghs[n], TPG); - } - storage_grad_hess_vector_[j] = dot(PG, storage_grad_hess_vector_[j]); - // grad_grad_grad_logdet(i,j)=storage_grad_hess_vector_[j]; - // grad_grad_psi(i,j)=storage_hess_vector_[j]; - // dpsi(i,j)=storage_grad_vector_[j]; - // psi(i,j)=storage_value_vector_[j]; - } - const std::complex eye(0.0, 1.0); - const std::complex meye(0.0, -1.0); - - // NumValenceOrbs appears to be the same as OrbitalSetSize in some code paths. - // In order to handle orbital rotation where OrbitalSetSize is not - // the same as the number of columns of psi, OrbitalSetSize was - // replaced with the columns of psi in other places. - // The minimum is used here just in case NumValenceOrbs might be smaller. - size_t loop_bound = NumValenceOrbs < psi.cols() ? NumValenceOrbs : psi.cols(); - for (int j = 0; j < loop_bound; j++) - { - std::complex u(storage_value_vector_[j]); - TinyVector, OHMMS_DIM> gradu(storage_grad_vector_[j]); - tmphs = storage_hess_vector_[j]; - tmpghs = storage_grad_hess_vector_[j]; - PosType k = kPoints[j]; - TinyVector ck; - for (int n = 0; n < OHMMS_DIM; n++) - ck[n] = k[n]; - double s, c; - double phase = -dot(r, k); - qmcplusplus::sincos(phase, &s, &c); - std::complex e_mikr(c, s); - psi(i, j) = e_mikr * u; - dpsi(i, j) = e_mikr * (-eye * u * ck + gradu); - grad_grad_psi(i, j) = - e_mikr * (tmphs - u * outerProduct(ck, ck) - eye * outerProduct(ck, gradu) - eye * outerProduct(gradu, ck)); - //Is this right? - storage_grad_hess_vector_[j] *= e_mikr; - for (unsigned a0(0); a0 < OHMMS_DIM; a0++) - for (unsigned a1(0); a1 < OHMMS_DIM; a1++) - for (unsigned a2(0); a2 < OHMMS_DIM; a2++) - storage_grad_hess_vector_[j][a0](a1, a2) += e_mikr * - (meye * (ck[a0] * tmphs(a1, a2) + ck[a1] * tmphs(a0, a2) + ck[a2] * tmphs(a0, a1)) - - (ck[a0] * ck[a1] * gradu[a2] + ck[a0] * ck[a2] * gradu[a1] + ck[a1] * ck[a2] * gradu[a0]) + - eye * ck[a0] * ck[a1] * ck[a2] * u); - grad_grad_grad_logdet(i, j) = storage_grad_hess_vector_[j]; - } - } -} - -#endif - -template -std::string EinsplineSetExtended::Type() -{ - return "EinsplineSetExtended"; -} - - -template -void EinsplineSetExtended::registerTimers() -{ - ValueTimer.reset(); - VGLTimer.reset(); - VGLMatTimer.reset(); - EinsplineTimer.reset(); -} - - -template -std::unique_ptr EinsplineSetExtended::makeClone() const -{ - auto clone = std::make_unique>(*this); - clone->registerTimers(); - for (int iat = 0; iat < clone->AtomicOrbitals.size(); iat++) - clone->AtomicOrbitals[iat].registerTimers(); - return clone; -} - -template class EinsplineSetExtended>; -template class EinsplineSetExtended; -} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/EinsplineSet.h b/src/QMCWaveFunctions/EinsplineSet.h deleted file mode 100644 index b4c65a0b8e..0000000000 --- a/src/QMCWaveFunctions/EinsplineSet.h +++ /dev/null @@ -1,359 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. -// -// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Ye Luo, yeluo@anl.gov, Argonne National Laboratory -// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory -// -// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -////////////////////////////////////////////////////////////////////////////////////// - - -#ifndef QMCPLUSPLUS_EINSPLINE_SET_H -#define QMCPLUSPLUS_EINSPLINE_SET_H - -#include "Configuration.h" -#include "QMCWaveFunctions/BasisSetBase.h" -#include "QMCWaveFunctions/SPOSet.h" -#include "QMCWaveFunctions/AtomicOrbital.h" -#include "Utilities/TimerManager.h" -#include "spline/einspline_engine.hpp" - -namespace qmcplusplus -{ -class EinsplineSetBuilder; - -class EinsplineSet : public SPOSet -{ - friend class EinsplineSetBuilder; - -public: - ////////////////////// - // Type definitions // - ////////////////////// - using UnitCellType = CrystalLattice; - - /////////// - // Flags // - /////////// - /// True if all Lattice is diagonal, i.e. 90 degree angles - bool Orthorhombic; - /// True if we are using localize orbitals - bool Localized; - /// True if we are tiling the primitive cell - bool Tiling; - - ////////////////////////// - // Lattice and geometry // - ////////////////////////// - TinyVector TileFactor; - Tensor TileMatrix; - UnitCellType SuperLattice, PrimLattice; - /// The "Twist" variables are in reduced coords, i.e. from 0 to1. - /// The "k" variables are in Cartesian coordinates. - PosType TwistVector, kVector; - /// This stores which "true" twist vector this clone is using. - /// "True" indicates the physical twist angle after untiling - int TwistNum; - /// metric tensor to handle generic unitcell - Tensor GGt; - - int NumValenceOrbs; - -public: - UnitCellType GetLattice(); - void resetSourceParticleSet(ParticleSet& ions); - void setOrbitalSetSize(int norbs) override; - inline std::string Type() { return "EinsplineSet"; } - EinsplineSet(const std::string& my_name) : SPOSet(my_name), TwistNum(0), NumValenceOrbs(0) {} - - virtual std::string getClassName() const override { return "EinsplineSet"; } -}; - -//////////////////////////////////////////////////////////////////// -// This is just a template trick to avoid template specialization // -// in EinsplineSetExtended. // -//////////////////////////////////////////////////////////////////// -template -struct MultiOrbitalTraits -{}; - -template<> -struct MultiOrbitalTraits -{ - using SplineType = multi_UBspline_2d_d; -}; - -template<> -struct MultiOrbitalTraits, 2> -{ - using SplineType = multi_UBspline_2d_z; -}; - -template<> -struct MultiOrbitalTraits -{ - using SplineType = multi_UBspline_2d_s; -}; - -template<> -struct MultiOrbitalTraits, 2> -{ - using SplineType = multi_UBspline_2d_c; -}; - -template<> -struct MultiOrbitalTraits -{ - using SplineType = multi_UBspline_3d_d; - using BCType = BCtype_d; - using DataType = double; -}; - -template<> -struct MultiOrbitalTraits, 3> -{ - using SplineType = multi_UBspline_3d_z; - using BCType = BCtype_z; - using DataType = std::complex; -}; - - -template<> -struct MultiOrbitalTraits -{ - using SplineType = multi_UBspline_3d_s; - using BCType = BCtype_s; - using DataType = float; -}; - -template<> -struct MultiOrbitalTraits, 3> -{ - using SplineType = multi_UBspline_3d_c; - using BCType = BCtype_c; - using DataType = std::complex; -}; - -//////////////////////////////////////////////////////////////////// -// Template class for evaluating multiple extended Bloch orbitals // -// quickly. Currently uses einspline library. // -//////////////////////////////////////////////////////////////////// -template -class EinsplineSetExtended : public EinsplineSet -{ - friend class EinsplineSetBuilder; - -protected: - ////////////////////// - // Type definitions // - ////////////////////// - //using UnitCellType = CrystalLattice; - using SplineType = typename MultiOrbitalTraits::SplineType; - using BCType = typename MultiOrbitalTraits::BCType; - - using StorageValueVector = typename OrbitalSetTraits::ValueVector; - using StorageGradVector = typename OrbitalSetTraits::GradVector; - using StorageHessVector = typename OrbitalSetTraits::HessVector; - using StorageGradHessVector = typename OrbitalSetTraits::GradHessVector; - using RealValueVector = Vector; - using ComplexValueVector = Vector>; - using RealGradVector = Vector>; - using ComplexGradVector = Vector, OHMMS_DIM>>; - using RealHessType = Tensor; - using ComplexHessType = Tensor, OHMMS_DIM>; - using RealHessVector = Vector; - using RealHessMatrix = Matrix; - using ComplexHessVector = Vector; - using ComplexHessMatrix = Matrix; - using RealValueMatrix = Matrix; - using ComplexValueMatrix = Matrix>; - using RealGradMatrix = Matrix>; - using ComplexGradMatrix = Matrix, OHMMS_DIM>>; - using RealGGGType = TinyVector; - using RealGGGVector = Vector; - using RealGGGMatrix = Matrix; - using ComplexGGGType = TinyVector; - using ComplexGGGVector = Vector; - using ComplexGGGMatrix = Matrix; - - ///////////////////////////// - /// Orbital storage object // - ///////////////////////////// - SplineType* MultiSpline; - - ////////////////////////////////////// - // Radial/Ylm orbitals around atoms // - ////////////////////////////////////// - std::vector> AtomicOrbitals; - - // First-order derivative w.r.t. the ion positions - std::vector> FirstOrderSplines; - // Temporary storage for Eispline calls - StorageValueVector storage_value_vector_, storage_lapl_vector_; - StorageGradVector storage_grad_vector_; - StorageHessVector storage_hess_vector_; - StorageGradHessVector storage_grad_hess_vector_; - - // True if we should unpack this orbital into two copies - std::vector MakeTwoCopies; - /** kpoints for each unique orbitals. - * Note: for historic reason, this sign is opposite to what was used in DFT when orbitals were generated. - * Changing the sign requires updating all the evaluation code. - */ - Vector> kPoints; - - /////////////////// - // Phase factors // - /////////////////// - Vector phase; - Vector> eikr; - void computePhaseFactors(const TinyVector& r); - // For running at half G-vectors with real orbitals; - // 0 if the twist is zero, 1 if the twist is G/2. - TinyVector HalfG; - - //////////// - // Timers // - //////////// - NewTimer &ValueTimer, &VGLTimer, &VGLMatTimer; - NewTimer& EinsplineTimer; - -public: - /** create MultiSpline - * @param xyz_g grid data - * @param xyz_bc boundary conditions - */ - template - void allocate(GT& xyz_g, BCT& xyz_bc, int nv) - { - SplineType* dummy = nullptr; - MultiSpline = einspline::create(dummy, xyz_g, xyz_bc, nv); - } - - inline void resizeStorage(int n, int nvals) - { - kPoints.resize(n); - MakeTwoCopies.resize(n); - storage_value_vector_.resize(n); - storage_lapl_vector_.resize(n); - storage_grad_vector_.resize(n); - storage_hess_vector_.resize(n); - storage_grad_hess_vector_.resize(n); - phase.resize(n); - eikr.resize(n); - NumValenceOrbs = nvals; - } - -#if !defined(QMC_COMPLEX) - // Real return values - void evaluateValue(const ParticleSet& P, int iat, RealValueVector& psi) override; - void evaluateVGL(const ParticleSet& P, - int iat, - RealValueVector& psi, - RealGradVector& dpsi, - RealValueVector& d2psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealValueMatrix& d2psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - RealValueMatrix& psi, - RealGradMatrix& dpsi, - RealHessMatrix& grad_grad_psi, - RealGGGMatrix& grad_grad_grad_logdet) override; - - // void evaluate (const ParticleSet& P, const PosType& r, std::vector &psi); - // This is the gradient of the orbitals w.r.t. the ion iat - void evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - RealGradMatrix& gradphi) override; - // Evaluate the gradient w.r.t. to ion iat of the gradient and - // laplacian of the orbitals w.r.t. the electrons - void evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - RealGradMatrix& dphi, - RealHessMatrix& dgrad_phi, - RealGradMatrix& dlaplphi) override; -#else - // Complex return values - void evaluateValue(const ParticleSet& P, int iat, ComplexValueVector& psi) override; - void evaluateVGL(const ParticleSet& P, - int iat, - ComplexValueVector& psi, - ComplexGradVector& dpsi, - ComplexValueVector& d2psi) override; - void evaluateVGH(const ParticleSet& P, - int iat, - ComplexValueVector& psi, - ComplexGradVector& dpsi, - ComplexHessVector& grad_grad_psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexValueMatrix& d2psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexHessMatrix& grad_grad_psi) override; - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ComplexValueMatrix& psi, - ComplexGradMatrix& dpsi, - ComplexHessMatrix& grad_grad_psi, - ComplexGGGMatrix& grad_grad_grad_logdet) override; -#endif - - void setOrbitalSetSize(int norbs) override; - std::string Type(); - - void registerTimers(); - PosType get_k(int orb) override { return kPoints[orb]; } - - virtual std::string getClassName() const override { return "EinsplineSetExtended"; } - bool hasIonDerivs() const override { return true; } - - std::unique_ptr makeClone() const override; - - EinsplineSetExtended(const std::string& my_name) - : EinsplineSet(my_name), - MultiSpline(NULL), - ValueTimer(createGlobalTimer("EinsplineSetExtended::ValueOnly")), - VGLTimer(createGlobalTimer("EinsplineSetExtended::VGL")), - VGLMatTimer(createGlobalTimer("EinsplineSetExtended::VGLMatrix")), - EinsplineTimer(createGlobalTimer("libeinspline")) - { - for (int i = 0; i < OHMMS_DIM; i++) - HalfG[i] = 0; - } -}; - -} // namespace qmcplusplus -#endif diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.h b/src/QMCWaveFunctions/EinsplineSetBuilder.h index edfaf9fc89..9e95677c48 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder.h +++ b/src/QMCWaveFunctions/EinsplineSetBuilder.h @@ -25,7 +25,6 @@ #include "QMCWaveFunctions/SPOSetBuilder.h" #include "QMCWaveFunctions/BandInfo.h" -#include "QMCWaveFunctions/AtomicOrbital.h" #include #include @@ -168,7 +167,6 @@ class EinsplineSetBuilder : public SPOSetBuilder TinyVector Version; std::string parameterGroup, ionsGroup, eigenstatesGroup; std::vector Occ; - bool HasCoreOrbs; bool ReadOrbitalInfo(bool skipChecks = false); bool ReadOrbitalInfo_ESHDF(bool skipChecks = false); void BroadcastOrbitalInfo(); @@ -197,7 +195,7 @@ class EinsplineSetBuilder : public SPOSetBuilder Tensor Lattice, RecipLattice, LatticeInv, SuperLattice, GGt; UnitCellType SuperCell, PrimCell, PrimCellInv; - int NumBands, NumElectrons, NumSpins, NumTwists, NumCoreStates; + int NumBands, NumElectrons, NumSpins, NumTwists; int MaxNumGvecs; double MeshFactor; RealType MatchingTol; @@ -226,7 +224,7 @@ class EinsplineSetBuilder : public SPOSetBuilder std::vector IncludeTwists, DistinctTwists; /// if false, splines are conceptually complex valued bool use_real_splines_; - int NumDistinctOrbitals, NumCoreOrbs, NumValenceOrbs; + int NumDistinctOrbitals; // This is true if the corresponding twist in DistinctTwists should // should be used to generate two distinct orbitals from the real and // imaginary parts. @@ -240,21 +238,9 @@ class EinsplineSetBuilder : public SPOSetBuilder void CopyBands(int numOrbs); - ///////////////////////////// - // Muffin-tin information // - ///////////////////////////// - int NumMuffinTins; - std::vector MT_APW_radii; - std::vector> MT_APW_rgrids; - std::vector MT_APW_lmax; - std::vector MT_APW_num_radial_points; - std::vector> MT_centers; - //////////////////////////////// // Atomic orbital information // //////////////////////////////// - std::vector>> AtomicOrbitals; - struct CenterInfo { std::vector lmax, spline_npoints, GroupID; @@ -281,8 +267,6 @@ class EinsplineSetBuilder : public SPOSetBuilder // This returns the path in the HDF5 file to the group for orbital // with twist ti and band bi std::string OrbitalPath(int ti, int bi); - std::string CoreStatePath(int ti, int bi); - std::string MuffinTinPath(int ti, int bi, int tin); ///////////////////////////////////////////////////////////// // Information to avoid storing the same orbitals twice in // @@ -298,9 +282,8 @@ class EinsplineSetBuilder : public SPOSetBuilder /** broadcast SortBands * @param N number of state * @param root true if it is the i/o node - * @return true, if core is found */ - bool bcastSortBands(int splin, int N, bool root); + void bcastSortBands(int splin, int N, bool root); /** a specific but clean code path in createSPOSetFromXML, for PBC, double, ESHDF * @param cur the current xml node diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index f35ffe73aa..054e27d949 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -48,12 +48,10 @@ EinsplineSetBuilder::EinsplineSetBuilder(ParticleSet& p, const PSetMap& psets, C NumElectrons(0), NumSpins(0), NumTwists(0), - NumCoreStates(0), MeshFactor(1.0), MeshSize(0, 0, 0), twist_num_(-1), TileFactor(1, 1, 1), - NumMuffinTins(0), LastSpinSet(-1), NumOrbitalsRead(-1), makeRotations(false) @@ -127,7 +125,6 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() if (myComm->size() == 1) return; int numIons = IonTypes.size(); - int numAtomicOrbitals = AtomicOrbitals.size(); int numDensityGvecs = TargetPtcl.DensityReducedGvecs.size(); PooledData abuffer; PooledData aibuffer; @@ -142,8 +139,6 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() aibuffer.add(NumSpins); //myComm->bcast(NumSpins); aibuffer.add(NumTwists); //myComm->bcast(NumTwists); aibuffer.add(numIons); //myComm->bcast(numIons); - aibuffer.add(NumMuffinTins); - aibuffer.add(numAtomicOrbitals); aibuffer.add(numDensityGvecs); aibuffer.add(HaveOrbDerivs); myComm->bcast(abuffer); @@ -163,26 +158,11 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() aibuffer.get(NumSpins); aibuffer.get(NumTwists); aibuffer.get(numIons); - aibuffer.get(NumMuffinTins); - aibuffer.get(numAtomicOrbitals); aibuffer.get(numDensityGvecs); aibuffer.get(HaveOrbDerivs); - MT_APW_radii.resize(NumMuffinTins); - MT_APW_lmax.resize(NumMuffinTins); - MT_APW_rgrids.resize(NumMuffinTins); - MT_APW_num_radial_points.resize(NumMuffinTins); - MT_centers.resize(NumMuffinTins); TargetPtcl.DensityReducedGvecs.resize(numDensityGvecs); TargetPtcl.Density_G.resize(numDensityGvecs); - AtomicOrbitals.resize(numAtomicOrbitals); } - std::vector rgrids_sizes(NumMuffinTins); - for (int tin = 0; tin < NumMuffinTins; tin++) - rgrids_sizes[tin] = MT_APW_rgrids[tin].size(); - myComm->bcast(rgrids_sizes); - if (myComm->rank()) - for (int tin = 0; tin < NumMuffinTins; tin++) - MT_APW_rgrids[tin].resize(rgrids_sizes[tin]); if (IonTypes.size() != numIons) { IonTypes.resize(numIons); @@ -199,27 +179,9 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() if (primcell_kpoints.size() != NumTwists) primcell_kpoints.resize(NumTwists); bbuffer.add(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); - bbuffer.add(MT_APW_radii.begin(), MT_APW_radii.end()); - bibuffer.add(MT_APW_lmax.begin(), MT_APW_lmax.end()); - bibuffer.add(MT_APW_num_radial_points.begin(), MT_APW_num_radial_points.end()); - bbuffer.add(&(MT_centers[0][0]), &(MT_centers[0][0]) + OHMMS_DIM * NumMuffinTins); - for (int i = 0; i < NumMuffinTins; i++) - bbuffer.add(MT_APW_rgrids[i].begin(), MT_APW_rgrids[i].end()); bibuffer.add(&(TargetPtcl.DensityReducedGvecs[0][0]), &(TargetPtcl.DensityReducedGvecs[0][0]) + numDensityGvecs * OHMMS_DIM); bbuffer.add(&(TargetPtcl.Density_G[0]), &(TargetPtcl.Density_G[0]) + numDensityGvecs); - for (int iat = 0; iat < numAtomicOrbitals; iat++) - { - AtomicOrbital>& orb = AtomicOrbitals[iat]; - bibuffer.add(orb.SplinePoints); - bibuffer.add(orb.PolyOrder); - bibuffer.add(orb.lMax); - bibuffer.add(orb.Numlm); - bbuffer.add(&orb.Pos[0], &orb.Pos[0] + OHMMS_DIM); - bbuffer.add(orb.CutoffRadius); - bbuffer.add(orb.SplineRadius); - bbuffer.add(orb.PolyRadius); - } myComm->bcast(bbuffer); myComm->bcast(bibuffer); if (myComm->rank()) @@ -230,27 +192,9 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() bibuffer.get(IonTypes[i]); bbuffer.get(&IonPos[0][0], &IonPos[0][0] + OHMMS_DIM * numIons); bbuffer.get(&primcell_kpoints[0][0], &primcell_kpoints[0][0] + OHMMS_DIM * NumTwists); - bbuffer.get(MT_APW_radii.begin(), MT_APW_radii.end()); - bibuffer.get(MT_APW_lmax.begin(), MT_APW_lmax.end()); - bibuffer.get(MT_APW_num_radial_points.begin(), MT_APW_num_radial_points.end()); - bbuffer.get(&(MT_centers[0][0]), &(MT_centers[0][0]) + OHMMS_DIM * NumMuffinTins); - for (int i = 0; i < NumMuffinTins; i++) - bbuffer.get(MT_APW_rgrids[i].begin(), MT_APW_rgrids[i].end()); bibuffer.get(&(TargetPtcl.DensityReducedGvecs[0][0]), &(TargetPtcl.DensityReducedGvecs[0][0]) + numDensityGvecs * OHMMS_DIM); bbuffer.get(&(TargetPtcl.Density_G[0]), &(TargetPtcl.Density_G[0]) + numDensityGvecs); - for (int iat = 0; iat < numAtomicOrbitals; iat++) - { - AtomicOrbital>& orb = AtomicOrbitals[iat]; - bibuffer.get(orb.SplinePoints); - bibuffer.get(orb.PolyOrder); - bibuffer.get(orb.lMax); - bibuffer.get(orb.Numlm); - bbuffer.get(&orb.Pos[0], &orb.Pos[0] + OHMMS_DIM); - bbuffer.get(orb.CutoffRadius); - bbuffer.get(orb.SplineRadius); - bbuffer.get(orb.PolyRadius); - } } //buffer to bcast hybrid representation atomic orbital info PooledData cbuffer; @@ -714,7 +658,6 @@ void EinsplineSetBuilder::OccupyBands(int spin, int sortBands, int numOrbs, bool for (int bi = 0; bi < NumBands; bi++) { BandInfo band; - band.IsCoreState = false; band.TwistIndex = tindex; band.BandIndex = bi; band.MakeTwoCopies = MakeTwoCopies[ti]; @@ -744,53 +687,31 @@ void EinsplineSetBuilder::OccupyBands(int spin, int sortBands, int numOrbs, bool SortBands.push_back(band); } } - // Now, read core states - for (int cs = 0; cs < NumCoreStates; cs++) - { - BandInfo band; - band.IsCoreState = true; - band.TwistIndex = tindex; - band.BandIndex = cs; - band.MakeTwoCopies = MakeTwoCopies[ti]; - H5File.read(band.Energy, CoreStatePath(ti, cs) + "eigenvalue"); - if (band.Energy > -1.0e100) - SortBands.push_back(band); - } } int orbIndex = 0; int numOrbs_counter = 0; - NumValenceOrbs = 0; - NumCoreOrbs = 0; while (numOrbs_counter < numOrbs) { if (SortBands[orbIndex].MakeTwoCopies) numOrbs_counter += 2; else numOrbs_counter++; - if (SortBands[orbIndex].IsCoreState) - NumCoreOrbs++; - else - NumValenceOrbs++; orbIndex++; } NumDistinctOrbitals = orbIndex; app_log() << "We will read " << NumDistinctOrbitals << " distinct orbitals.\n"; - app_log() << "There are " << NumCoreOrbs << " core states and " << NumValenceOrbs << " valence states.\n"; } -bool EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) +void EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) { std::vector& SortBands(*FullBands[spin]); - TinyVector nbands(int(SortBands.size()), n, NumValenceOrbs, NumCoreOrbs); + TinyVector nbands(int(SortBands.size()), n); mpi::bcast(*myComm, nbands); //buffer to serialize BandInfo - PooledData misc(nbands[0] * 5); - bool isCore = false; + PooledData misc(nbands[0] * 4); n = NumDistinctOrbitals = nbands[1]; - NumValenceOrbs = nbands[2]; - NumCoreOrbs = nbands[3]; if (root) { @@ -801,9 +722,6 @@ bool EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) misc.put(SortBands[i].BandIndex); misc.put(SortBands[i].Energy); misc.put(SortBands[i].MakeTwoCopies); - misc.put(SortBands[i].IsCoreState); - - isCore |= SortBands[i].IsCoreState; } for (int i = n; i < SortBands.size(); ++i) @@ -812,7 +730,6 @@ bool EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) misc.put(SortBands[i].BandIndex); misc.put(SortBands[i].Energy); misc.put(SortBands[i].MakeTwoCopies); - misc.put(SortBands[i].IsCoreState); } } myComm->bcast(misc); @@ -827,9 +744,6 @@ bool EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) misc.get(SortBands[i].BandIndex); misc.get(SortBands[i].Energy); misc.get(SortBands[i].MakeTwoCopies); - misc.get(SortBands[i].IsCoreState); - - isCore |= SortBands[i].IsCoreState; } for (int i = n; i < SortBands.size(); ++i) { @@ -837,10 +751,8 @@ bool EinsplineSetBuilder::bcastSortBands(int spin, int n, bool root) misc.get(SortBands[i].BandIndex); misc.get(SortBands[i].Energy); misc.get(SortBands[i].MakeTwoCopies); - misc.get(SortBands[i].IsCoreState); } } - return isCore; } } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp index 32f547a50b..5f855d6257 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp @@ -76,20 +76,15 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) for (int j = 0; j < 3; j++) LatticeInv(i, j) = RecipLattice(i, j) / (2.0 * M_PI); int have_dpsi = false; - int NumAtomicOrbitals = 0; - NumCoreStates = NumMuffinTins = NumTwists = NumSpins = NumBands = NumAtomicOrbitals = 0; + NumTwists = NumSpins = NumBands = 0; NumElectrons = TargetPtcl.getTotalNum(); H5File.read(NumBands, "/electrons/kpoint_0/spin_0/number_of_states"); - H5File.readEntry(NumCoreStates, "/electrons/kpoint_0/spin_0/number_of_core_states"); H5File.readEntry(NumSpins, "/electrons/number_of_spins"); H5File.read(NumTwists, "/electrons/number_of_kpoints"); - H5File.readEntry(NumMuffinTins, "/muffin_tins/number_of_tins"); H5File.readEntry(have_dpsi, "/electrons/have_dpsi"); - H5File.readEntry(NumAtomicOrbitals, "/electrons/number_of_atomic_orbitals"); HaveOrbDerivs = have_dpsi; app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists - << ", muffin tins=" << NumMuffinTins << ", core states=" << NumCoreStates << std::endl; - app_log() << "atomic orbital=" << NumAtomicOrbitals << std::endl; + << std::endl; if (TileFactor[0] != 1 || TileFactor[1] != 1 || TileFactor[2] != 1) app_log() << " Using a " << TileFactor[0] << "x" << TileFactor[1] << "x" << TileFactor[2] << " tiling factor.\n"; ////////////////////////////////// @@ -229,38 +224,6 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) AtomicCentersInfo.lmax[center_idx] = source_species(lmax_ind, my_GroupID); } } - ///////////////////////////////////// - // Read atomic orbital information // - ///////////////////////////////////// - AtomicOrbitals.resize(NumAtomicOrbitals); - for (int iat = 0; iat < NumAtomicOrbitals; iat++) - { - AtomicOrbital>& orb = AtomicOrbitals[iat]; - int lmax, polynomial_order, spline_points; - RealType cutoff_radius, polynomial_radius, spline_radius; - PosType position; - double cutoff_radius_DP, polynomial_radius_DP, spline_radius_DP; - TinyVector position_DP; - std::ostringstream groupstream; - groupstream << "/electrons/atomic_orbital_" << iat << "/"; - std::string groupname = groupstream.str(); - H5File.read(lmax, groupname + "lmax"); - H5File.read(polynomial_order, groupname + "polynomial_order"); - H5File.read(spline_points, groupname + "spline_points"); - H5File.read(cutoff_radius_DP, groupname + "cutoff_radius"); - H5File.read(polynomial_radius_DP, groupname + "polynomial_radius"); - H5File.read(spline_radius_DP, groupname + "spline_radius"); - H5File.read(position_DP, groupname + "position"); - cutoff_radius = cutoff_radius_DP; - polynomial_radius = polynomial_radius_DP; - spline_radius = spline_radius_DP; - position = position_DP; - orb.set_pos(position); - orb.set_lmax(lmax); - orb.set_cutoff(cutoff_radius); - orb.set_spline(spline_radius, spline_points); - orb.set_polynomial(polynomial_radius, polynomial_order); - } /////////////////////////// // Read the twist angles // /////////////////////////// @@ -398,7 +361,6 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs for (int bi = 0; bi < NumBands; bi++) { BandInfo band; - band.IsCoreState = false; band.TwistIndex = tindex; band.BandIndex = bi; band.MakeTwoCopies = MakeTwoCopies[ti]; @@ -410,22 +372,6 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs else maxOrbs++; } - // Now, read core states - for (int cs = 0; cs < NumCoreStates; cs++) - { - BandInfo band; - band.IsCoreState = true; - band.TwistIndex = tindex; - band.BandIndex = cs; - band.MakeTwoCopies = MakeTwoCopies[ti]; - H5File.read(band.Energy, CoreStatePath(ti, cs) + "eigenvalue"); - if (band.Energy > -1.0e100) - SortBands.push_back(band); - if (MakeTwoCopies[ti]) - maxOrbs += 2; - else - maxOrbs++; - } } app_log() << SortBands.size() << " complex-valued orbitals supplied by h5 can be expanded up to " << maxOrbs @@ -594,23 +540,16 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs //} int orbIndex = 0; int numOrbs_counter = 0; - NumValenceOrbs = 0; - NumCoreOrbs = 0; while (numOrbs_counter < numOrbs) { if (SortBands[orbIndex].MakeTwoCopies) numOrbs_counter += 2; else numOrbs_counter++; - if (SortBands[orbIndex].IsCoreState) - NumCoreOrbs++; - else - NumValenceOrbs++; orbIndex++; } NumDistinctOrbitals = orbIndex; app_log() << "We will read " << NumDistinctOrbitals << " distinct complex-valued orbitals from h5.\n"; - app_log() << "There are " << NumCoreOrbs << " core states and " << NumValenceOrbs << " valence states.\n"; } } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp index c4d316445a..0ba6be3cb0 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp @@ -97,39 +97,15 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) LatticeInv(i, j) = RecipLattice(i, j) / (2.0 * M_PI); - NumCoreStates = NumMuffinTins = 0; H5File.read(NumBands, parameterGroup + "/num_bands"); - H5File.read(NumCoreStates, parameterGroup + "/num_core_states"); H5File.read(NumElectrons, parameterGroup + "/num_electrons"); H5File.read(NumSpins, parameterGroup + "/num_spins"); H5File.read(NumTwists, parameterGroup + "/num_twists"); - H5File.read(NumMuffinTins, parameterGroup + "/muffin_tins/num_tins"); app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists - << ", muffin tins=" << NumMuffinTins << std::endl; + << std::endl; if (TileFactor[0] != 1 || TileFactor[1] != 1 || TileFactor[2] != 1) app_log() << " Using a " << TileFactor[0] << "x" << TileFactor[1] << "x" << TileFactor[2] << " tiling factor.\n"; - ///////////////////////////////// - // Read muffin tin information // - ///////////////////////////////// - MT_APW_radii.resize(NumMuffinTins); - MT_APW_rgrids.resize(NumMuffinTins); - MT_APW_lmax.resize(NumMuffinTins); - MT_APW_num_radial_points.resize(NumMuffinTins); - MT_centers.resize(NumMuffinTins); - for (int tin = 0; tin < NumMuffinTins; tin++) - { - std::ostringstream MTstream; - if (NumMuffinTins > 1) - MTstream << parameterGroup << "/muffin_tins/muffin_tin_" << tin; - else - MTstream << parameterGroup << "/muffin_tins/muffin_tin"; - std::string MTgroup = MTstream.str(); - H5File.read(MT_APW_lmax[tin], MTgroup + "/lmax"); - H5File.read(MT_APW_num_radial_points[tin], MTgroup + "/num_radial_points"); - H5File.read(MT_APW_radii[tin], MTgroup + "/radius"); - H5File.read(MT_centers[tin], MTgroup + "/center"); - H5File.read(MT_APW_rgrids[tin], MTgroup + "/r"); - } + ////////////////////////////////// // Read ion types and locations // ////////////////////////////////// @@ -188,48 +164,4 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) return true; } - -std::string EinsplineSetBuilder::OrbitalPath(int ti, int bi) -{ - std::string eigenstatesGroup; - if (Version[0] == 0 && Version[1] == 11) - eigenstatesGroup = "/eigenstates_3"; - else if (Version[0] == 0 && Version[1] == 20) - eigenstatesGroup = "/eigenstates"; - std::ostringstream groupPath; - if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1) - groupPath << eigenstatesGroup << "/twist_" << ti << "/band_" << bi << "/"; - else if (NumBands > 1) - groupPath << eigenstatesGroup << "/twist/band_" << bi << "/"; - else - groupPath << eigenstatesGroup << "/twist/band/"; - return groupPath.str(); -} - -std::string EinsplineSetBuilder::CoreStatePath(int ti, int cs) -{ - std::string eigenstatesGroup; - if (Version[0] == 0 && Version[1] == 11) - eigenstatesGroup = "/eigenstates_3"; - else if (Version[0] == 0 && Version[1] == 20) - eigenstatesGroup = "/eigenstates"; - std::ostringstream groupPath; - if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1) - groupPath << eigenstatesGroup << "/twist_" << ti << "/core_state_" << cs << "/"; - else if (NumBands > 1) - groupPath << eigenstatesGroup << "/twist/core_state_" << cs << "/"; - else - groupPath << eigenstatesGroup << "/twist/core_state/"; - return groupPath.str(); -} - -std::string EinsplineSetBuilder::MuffinTinPath(int ti, int bi, int tin) -{ - std::ostringstream groupPath; - if (NumMuffinTins > 0) - groupPath << OrbitalPath(ti, bi) << "muffin_tin_" << tin << "/"; - else - groupPath << OrbitalPath(ti, bi) << "muffin_tin/"; - return groupPath.str(); -} } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp index 223e689f58..e610b6bac8 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp @@ -29,9 +29,6 @@ #include #include "Utilities/ProgressReportEngine.h" #include "QMCWaveFunctions/einspline_helper.hpp" -#if !defined(MIXED_PRECISION) -#include "QMCWaveFunctions/EinsplineSet.h" -#endif #include "QMCWaveFunctions/BsplineFactory/BsplineReaderBase.h" #include "QMCWaveFunctions/BsplineFactory/BsplineSet.h" #include "QMCWaveFunctions/BsplineFactory/createBsplineReader.h" @@ -104,8 +101,6 @@ void EinsplineSetBuilder::set_metadata(int numOrbs, PrimCell.set(Lattice); SuperCell.set(SuperLattice); GGt = dot(transpose(PrimCell.G), PrimCell.G); - for (int iat = 0; iat < AtomicOrbitals.size(); iat++) - AtomicOrbitals[iat].Lattice.set(Lattice); // Now, analyze the k-point mesh to figure out the what k-points are needed AnalyzeTwists2(twist_num_inp, twist_inp); @@ -149,7 +144,6 @@ std::unique_ptr EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur) a.add(GPUsharing, "gpusharing"); // split spline across GPUs visible per rank a.add(spo_prec, "precision"); a.add(truncate, "truncate"); - a.add(use_einspline_set_extended, "use_old_spline"); a.add(myName, "tag"); a.add(skip_checks, "skip_checks"); @@ -270,9 +264,6 @@ std::unique_ptr EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur) kid = kid->next; } - if (has_backflow && use_einspline_set_extended == "yes" && use_real_splines_) - myComm->barrier_and_abort("backflow optimization is broken with use_real_splines_"); - ////////////////////////////////// // Create the OrbitalSet object ////////////////////////////////// @@ -316,104 +307,11 @@ std::unique_ptr EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur) MixedSplineReader->setCommon(XMLRoot); // temporary disable the following function call, Ye Luo // RotateBands_ESHDF(spinSet, dynamic_cast >*>(OrbitalSet)); - HasCoreOrbs = bcastSortBands(spinSet, NumDistinctOrbitals, myComm->rank() == 0); + bcastSortBands(spinSet, NumDistinctOrbitals, myComm->rank() == 0); auto OrbitalSet = MixedSplineReader->create_spline_set(spinSet, spo_cur); if (!OrbitalSet) myComm->barrier_and_abort("Failed to create SPOSet*"); -#if defined(MIXED_PRECISION) - if (use_einspline_set_extended == "yes") - myComm->barrier_and_abort("Option use_old_spline is not supported by the mixed precision build!"); -#else - if (use_einspline_set_extended == "yes") - { - std::unique_ptr new_OrbitalSet; - if (use_real_splines_) - { - auto temp_OrbitalSet = std::make_unique>(spo_object_name); - temp_OrbitalSet->MultiSpline = MixedSplineReader->export_MultiSplineDouble().release(); - temp_OrbitalSet->MultiSpline->num_splines = NumDistinctOrbitals; - temp_OrbitalSet->resizeStorage(NumDistinctOrbitals, NumValenceOrbs); - //set the flags for anti periodic boundary conditions - temp_OrbitalSet->HalfG = dynamic_cast(*OrbitalSet).getHalfG(); - new_OrbitalSet = std::move(temp_OrbitalSet); - } - else - { - auto temp_OrbitalSet = std::make_unique>>(spo_object_name); - temp_OrbitalSet->MultiSpline = MixedSplineReader->export_MultiSplineComplexDouble().release(); - temp_OrbitalSet->MultiSpline->num_splines = NumDistinctOrbitals; - temp_OrbitalSet->resizeStorage(NumDistinctOrbitals, NumValenceOrbs); - for (int iorb = 0, num = 0; iorb < NumDistinctOrbitals; iorb++) - { - int ti = (*FullBands[spinSet])[iorb].TwistIndex; - temp_OrbitalSet->kPoints[iorb] = PrimCell.k_cart(-primcell_kpoints[ti]); - temp_OrbitalSet->MakeTwoCopies[iorb] = (num < (numOrbs - 1)) && (*FullBands[spinSet])[iorb].MakeTwoCopies; - num += temp_OrbitalSet->MakeTwoCopies[iorb] ? 2 : 1; - } - new_OrbitalSet = std::move(temp_OrbitalSet); - } - //set the internal parameters - setTiling(new_OrbitalSet.get(), numOrbs); - OrbitalSet = std::move(new_OrbitalSet); - } -#endif app_log() << "Time spent in creating B-spline SPOs " << mytimer.elapsed() << "sec" << std::endl; -#ifdef Ye_debug -#ifndef QMC_COMPLEX - if (myComm->rank() == 0 && OrbitalSet->MuffinTins.size() > 0) - { - FILE* fout = fopen("TestMuffins.dat", "w"); - Vector phi(numOrbs), lapl(numOrbs); - Vector grad(numOrbs); - ParticleSet P; - P.R.resize(6); - for (int i = 0; i < P.R.size(); i++) - P.R[i] = PosType(0.0, 0.0, 0.0); - PosType N = 0.25 * PrimCell.a(0) + 0.25 * PrimCell.a(1) + 0.25 * PrimCell.a(2); - for (double x = -1.0; x <= 1.0; x += 0.0000500113412) - { - // for (double x=-0.003; x<=0.003; x+=0.0000011329343481381) { - P.R[0] = x * (PrimCell.a(0) + 0.914 * PrimCell.a(1) + 0.781413 * PrimCell.a(2)); - double r = std::sqrt(dot(P.R[0], P.R[0])); - double rN = std::sqrt(dot(P.R[0] - N, P.R[0] - N)); - OrbitalSet->evaluate(P, 0, phi, grad, lapl); - // OrbitalSet->evaluate(P, 0, phi); - fprintf(fout, "%1.12e ", r * x / std::abs(x)); - for (int j = 0; j < numOrbs; j++) - { - double gmag = std::sqrt(dot(grad[j], grad[j])); - fprintf(fout, "%16.12e ", - /*phi[j]*phi[j]**/ (-5.0 / r - 0.5 * lapl[j] / phi[j])); - // double E = -5.0/r -0.5*lapl[j]/phi[j]; - fprintf(fout, "%16.12e ", phi[j]); - fprintf(fout, "%16.12e ", gmag); - } - fprintf(fout, "\n"); - } - fclose(fout); - } -#endif -#endif - //if (sourceName.size() && (ParticleSets.find(sourceName) == ParticleSets.end())) - //{ - // app_log() << " EinsplineSetBuilder creates a ParticleSet " << sourceName << std::endl; - // ParticleSet* ions=new ParticleSet; - // ions->Lattice=TargetPtcl.Lattice; - // ESHDFIonsParser ap(*ions,H5FileID,myComm); - // ap.put(XMLRoot); - // ap.expand(TileMatrix); - // ions->setName(sourceName); - // ParticleSets[sourceName]=ions; - // //overwrite the lattice and assign random - // if(TargetPtcl.Lattice.SuperCellEnum) - // { - // TargetPtcl.Lattice=ions->Lattice; - // makeUniformRandom(TargetPtcl.R); - // TargetPtcl.R.setUnit(PosUnit::LatticeUnit); - // TargetPtcl.convert2Cart(TargetPtcl.R); - // TargetPtcl.createSK(); - // } - //} OrbitalSet->finalizeConstruction(); SPOSetMap[aset] = OrbitalSet.get(); return OrbitalSet; diff --git a/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp b/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp index 85a2bf4509..c028ccc630 100644 --- a/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp +++ b/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp @@ -199,12 +199,12 @@ std::unique_ptr EinsplineSpinorSetBuilder::createSPOSetFromXML(xmlNodePt MixedSplineReader->setRotate(false); //Make the up spin set. - HasCoreOrbs = bcastSortBands(spinSet, NumDistinctOrbitals, myComm->rank() == 0); + bcastSortBands(spinSet, NumDistinctOrbitals, myComm->rank() == 0); auto bspline_zd_u = MixedSplineReader->create_spline_set(spinSet, spo_cur); //Make the down spin set. OccupyBands(spinSet2, sortBands, numOrbs, skipChecks); - HasCoreOrbs = bcastSortBands(spinSet2, NumDistinctOrbitals, myComm->rank() == 0); + bcastSortBands(spinSet2, NumDistinctOrbitals, myComm->rank() == 0); auto bspline_zd_d = MixedSplineReader->create_spline_set(spinSet2, spo_cur); //register with spin set and we're off to the races. From b1de0bc997eaf7f91ef017af32aea29da068a825 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 23:06:28 -0500 Subject: [PATCH 5/9] Remove TileFactor --- src/QMCWaveFunctions/EinsplineSetBuilder.h | 28 +++---------------- .../EinsplineSetBuilderCommon.cpp | 3 +- .../EinsplineSetBuilderESHDF.fft.cpp | 2 -- .../EinsplineSetBuilderOld.cpp | 2 -- .../EinsplineSetBuilder_createSPOs.cpp | 8 +++--- 5 files changed, 9 insertions(+), 34 deletions(-) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.h b/src/QMCWaveFunctions/EinsplineSetBuilder.h index 9e95677c48..f7794dbba4 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder.h +++ b/src/QMCWaveFunctions/EinsplineSetBuilder.h @@ -177,22 +177,6 @@ class EinsplineSetBuilder : public SPOSetBuilder */ bool ReadGvectors_ESHDF(); - /** set tiling properties of oset - * @param oset spline-orbital engine to be initialized - * @param numOrbs number of orbitals that belong to oset - */ - template - inline void setTiling(SPE* oset, int numOrbs) - { - oset->TileFactor = TileFactor; - oset->Tiling = (TileFactor[0] * TileFactor[1] * TileFactor[2] != 1); - oset->PrimLattice.set(Lattice); - oset->SuperLattice.set(SuperLattice); - oset->GGt = GGt; - oset->setOrbitalSetSize(numOrbs); - } - - Tensor Lattice, RecipLattice, LatticeInv, SuperLattice, GGt; UnitCellType SuperCell, PrimCell, PrimCellInv; int NumBands, NumElectrons, NumSpins, NumTwists; @@ -214,12 +198,9 @@ class EinsplineSetBuilder : public SPOSetBuilder int twist_num_; // primitive cell k-points from DFT calculations std::vector> primcell_kpoints; - - TinyVector TileFactor; + // primitive cell to supercell tiling matrix Tensor TileMatrix; - TinyVector TwistMesh; - // This vector stores which twist indices will be used by this - // clone + // This vector stores which twist indices will be used by this clone std::vector> UseTwists; std::vector IncludeTwists, DistinctTwists; /// if false, splines are conceptually complex valued @@ -229,15 +210,14 @@ class EinsplineSetBuilder : public SPOSetBuilder // should be used to generate two distinct orbitals from the real and // imaginary parts. std::vector MakeTwoCopies; - inline bool TwistPair(PosType a, PosType b); // This maps a 3-integer twist index into the twist number in the file std::map, int, Int3less> TwistMap; + + bool TwistPair(PosType a, PosType b) const; void TileIons(); void OccupyBands(int spin, int sortBands, int numOrbs, bool skipChecks = false); void OccupyBands_ESHDF(int spin, int sortBands, int numOrbs); - void CopyBands(int numOrbs); - //////////////////////////////// // Atomic orbital information // //////////////////////////////// diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index 054e27d949..101b587317 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -51,7 +51,6 @@ EinsplineSetBuilder::EinsplineSetBuilder(ParticleSet& p, const PSetMap& psets, C MeshFactor(1.0), MeshSize(0, 0, 0), twist_num_(-1), - TileFactor(1, 1, 1), LastSpinSet(-1), NumOrbitalsRead(-1), makeRotations(false) @@ -319,7 +318,7 @@ void EinsplineSetBuilder::TileIons() } -bool EinsplineSetBuilder::TwistPair(PosType a, PosType b) +bool EinsplineSetBuilder::TwistPair(PosType a, PosType b) const { bool pair = true; for (int n = 0; n < OHMMS_DIM; n++) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp index 5f855d6257..0c8d9095a7 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp @@ -85,8 +85,6 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) HaveOrbDerivs = have_dpsi; app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists << std::endl; - if (TileFactor[0] != 1 || TileFactor[1] != 1 || TileFactor[2] != 1) - app_log() << " Using a " << TileFactor[0] << "x" << TileFactor[1] << "x" << TileFactor[2] << " tiling factor.\n"; ////////////////////////////////// // Read ion types and locations // ////////////////////////////////// diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp index 0ba6be3cb0..7d10e81eab 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp @@ -103,8 +103,6 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) H5File.read(NumTwists, parameterGroup + "/num_twists"); app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists << std::endl; - if (TileFactor[0] != 1 || TileFactor[1] != 1 || TileFactor[2] != 1) - app_log() << " Using a " << TileFactor[0] << "x" << TileFactor[1] << "x" << TileFactor[2] << " tiling factor.\n"; ////////////////////////////////// // Read ion types and locations // diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp index e610b6bac8..9c02d343a9 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilder_createSPOs.cpp @@ -58,12 +58,11 @@ void EinsplineSetBuilder::set_metadata(int numOrbs, for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) matrixNotSet = matrixNotSet && (TileMatrix(i, j) == 0); - // then set the matrix to what may have been specified in the - // tiling vector + // then set the matrix to identity. if (matrixNotSet) for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) - TileMatrix(i, j) = (i == j) ? TileFactor[i] : 0; + TileMatrix(i, j) = (i == j) ? 1 : 0; if (myComm->rank() == 0) { std::array buff; @@ -130,9 +129,10 @@ std::unique_ptr EinsplineSetBuilder::createSPOSetFromXML(xmlNodePtr cur) ScopedTimer spo_timer_scope(createGlobalTimer("einspline::CreateSPOSetFromXML", timer_level_medium)); { + TinyVector TileFactor_do_not_use; OhmmsAttributeSet a; a.add(H5FileName, "href"); - a.add(TileFactor, "tile"); + a.add(TileFactor_do_not_use, "tile", {}, TagStatus::DELETED); a.add(sortBands, "sort"); a.add(TileMatrix, "tilematrix"); a.add(twist_num_inp, "twistnum"); From f4361925c5db3879e380c4ad427e188c7b704264 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 23:34:53 -0500 Subject: [PATCH 6/9] cut down EinsplineSetBuilder::ReadOrbitalInfo --- .../EinsplineSetBuilderOld.cpp | 121 +----------------- 1 file changed, 6 insertions(+), 115 deletions(-) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp index 7d10e81eab..2e84a56383 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderOld.cpp @@ -33,10 +33,10 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) { if (!H5File.open(H5FileName, H5F_ACC_RDONLY)) { - app_error() << "Could not open HDF5 file \"" << H5FileName - << "\" in EinsplineSetBuilder::ReadOrbitalInfo. Aborting.\n"; - APP_ABORT("EinsplineSetBuilder::ReadOrbitalInfo"); + app_error() << "Could not open HDF5 file \"" << H5FileName << "\" in EinsplineSetBuilder::ReadOrbitalInfo.\n"; + return false; } + // Read format std::string format; H5File.read(format, "/format"); @@ -47,119 +47,10 @@ bool EinsplineSetBuilder::ReadOrbitalInfo(bool skipChecks) Format = ESHDF; return ReadOrbitalInfo_ESHDF(skipChecks); } - ////////////////////////////////////////////////// - // Read basic parameters from the orbital file. // - ////////////////////////////////////////////////// - // Check the version - if (Version[0] == 0 && Version[1] == 11) - { - parameterGroup = "/parameters_0"; - ionsGroup = "/ions_2"; - eigenstatesGroup = "/eigenstates_3"; - } - else if (Version[0] == 0 && Version[1] == 20) - { - parameterGroup = "/parameters"; - ionsGroup = "/ions"; - eigenstatesGroup = "/eigenstates"; - } - else - { - std::ostringstream o; - o << "Unknown HDF5 orbital file version " << Version[0] << "." << Version[1] << "." << Version[2] << "\n"; - APP_ABORT(o.str()); - } - H5File.read(Lattice, parameterGroup + "/lattice"); - H5File.read(RecipLattice, parameterGroup + "/reciprocal_lattice"); - SuperLattice = dot(TileMatrix, Lattice); - std::array buff; - int length = std::snprintf(buff.data(), buff.size(), - " Lattice = \n [ %8.5f %8.5f %8.5f\n" - " %8.5f %8.5f %8.5f\n" - " %8.5f %8.5f %8.5f ]\n", - Lattice(0, 0), Lattice(0, 1), Lattice(0, 2), Lattice(1, 0), Lattice(1, 1), Lattice(1, 2), - Lattice(2, 0), Lattice(2, 1), Lattice(2, 2)); - if (length < 0) - throw std::runtime_error("Error generating Lattice string"); - app_log() << std::string_view(buff.data(), length); - length = - std::snprintf(buff.data(), buff.size(), - " SuperLattice = \n [ %13.12f %13.12f %13.12f\n" - " %13.12f %13.12f %13.12f\n" - " %13.12f %13.12f %13.12f ]\n", - SuperLattice(0, 0), SuperLattice(0, 1), SuperLattice(0, 2), SuperLattice(1, 0), SuperLattice(1, 1), - SuperLattice(1, 2), SuperLattice(2, 0), SuperLattice(2, 1), SuperLattice(2, 2)); - if (length < 0) - throw std::runtime_error("Error generating SuperLattice string"); - if (!CheckLattice()) - APP_ABORT("CheckLattice failed"); - app_log() << std::string_view(buff.data(), length); - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) - LatticeInv(i, j) = RecipLattice(i, j) / (2.0 * M_PI); - H5File.read(NumBands, parameterGroup + "/num_bands"); - H5File.read(NumElectrons, parameterGroup + "/num_electrons"); - H5File.read(NumSpins, parameterGroup + "/num_spins"); - H5File.read(NumTwists, parameterGroup + "/num_twists"); - app_log() << "bands=" << NumBands << ", elecs=" << NumElectrons << ", spins=" << NumSpins << ", twists=" << NumTwists - << std::endl; - ////////////////////////////////// - // Read ion types and locations // - ////////////////////////////////// - H5File.read(IonTypes, ionsGroup + "/atom_types"); - H5File.read(IonPos, ionsGroup + "/pos"); - /////////////////////////// - // Read the twist angles // - /////////////////////////// - primcell_kpoints.resize(NumTwists); - for (int ti = 0; ti < NumTwists; ti++) - { - std::ostringstream path; - if ((Version[0] == 0 && Version[1] == 11) || NumTwists > 1) - path << eigenstatesGroup << "/twist_" << ti << "/twist_angle"; - else - path << eigenstatesGroup << "/twist/twist_angle"; - TinyVector primcell_kpoints_DP; - H5File.read(primcell_kpoints_DP, path.str()); - primcell_kpoints[ti] = primcell_kpoints_DP; - int length = std::snprintf(buff.data(), buff.size(), " Found twist angle (%6.3f, %6.3f, %6.3f)\n", - primcell_kpoints[ti][0], primcell_kpoints[ti][1], primcell_kpoints[ti][2]); - if (length < 0) - throw std::runtime_error("Error converting twist angle to string"); - app_log() << std::string_view(buff.data(), length); - } - ////////////////////////////////////////////////////////// - // If the density has not been set in TargetPtcl, and // - // the density is available, read it in and save it // - // in TargetPtcl. // - ////////////////////////////////////////////////////////// - if (TargetPtcl.Density_G.empty()) - { - Array Density_r_DP; - H5File.read(TargetPtcl.DensityReducedGvecs, "/density/reduced_gvecs"); - H5File.read(Density_r_DP, "/density/rho_r"); - TargetPtcl.Density_r = Density_r_DP; - int numG = TargetPtcl.DensityReducedGvecs.size(); - // Convert primitive G-vectors to supercell G-vectors - for (int iG = 0; iG < numG; iG++) - TargetPtcl.DensityReducedGvecs[iG] = dot(TileMatrix, TargetPtcl.DensityReducedGvecs[iG]); - app_log() << " Read " << numG << " density G-vectors.\n"; - if (TargetPtcl.DensityReducedGvecs.size()) - { - app_log() << " EinsplineSetBuilder found density in the HDF5 file.\n"; - std::vector> Density_G_DP; - H5File.read(Density_G_DP, "/density/rho_G"); - TargetPtcl.Density_G.assign(Density_G_DP.begin(), Density_G_DP.end()); - if (!TargetPtcl.Density_G.size()) - { - app_error() << " Density reduced G-vectors defined, but not the" - << " density.\n"; - abort(); - } - } - } - return true; + app_error() + << "EinsplineSetBuilder::ReadOrbitalInfo too old h5 file which is not in ESHDF format! Regenerate the h5 file"; + return false; } } // namespace qmcplusplus From 8de1b42b4a646a61a3d04983ae7d8b829938d548 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Thu, 27 Jul 2023 23:35:41 -0500 Subject: [PATCH 7/9] Formatting. --- src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp | 4 ++-- src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp | 7 ++++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp index 101b587317..04f2492455 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderCommon.cpp @@ -123,8 +123,8 @@ void EinsplineSetBuilder::BroadcastOrbitalInfo() { if (myComm->size() == 1) return; - int numIons = IonTypes.size(); - int numDensityGvecs = TargetPtcl.DensityReducedGvecs.size(); + int numIons = IonTypes.size(); + int numDensityGvecs = TargetPtcl.DensityReducedGvecs.size(); PooledData abuffer; PooledData aibuffer; aibuffer.add(Version.begin(), Version.end()); //myComm->bcast(Version); diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp index 0c8d9095a7..82c4507df6 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderESHDF.fft.cpp @@ -75,9 +75,9 @@ bool EinsplineSetBuilder::ReadOrbitalInfo_ESHDF(bool skipChecks) for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) LatticeInv(i, j) = RecipLattice(i, j) / (2.0 * M_PI); - int have_dpsi = false; + int have_dpsi = false; NumTwists = NumSpins = NumBands = 0; - NumElectrons = TargetPtcl.getTotalNum(); + NumElectrons = TargetPtcl.getTotalNum(); H5File.read(NumBands, "/electrons/kpoint_0/spin_0/number_of_states"); H5File.readEntry(NumSpins, "/electrons/number_of_spins"); H5File.read(NumTwists, "/electrons/number_of_kpoints"); @@ -411,7 +411,8 @@ void EinsplineSetBuilder::OccupyBands_ESHDF(int spin, int sortBands, int numOrbs } if (occ_format == "energy") { - app_log() << " Occupying bands based on energy in mode " << (Occ.size() > 0? "\"excited\"" : "\"ground\"") << std::endl; + app_log() << " Occupying bands based on energy in mode " << (Occ.size() > 0 ? "\"excited\"" : "\"ground\"") + << std::endl; // To get the occupations right. std::vector Removed(0, 0); std::vector Added(0, 0); From 6910bd7c92f953b8a665732fe28424e8e4f470a8 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Fri, 28 Jul 2023 13:38:21 -0500 Subject: [PATCH 8/9] Fix complex build. --- src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp b/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp index c028ccc630..a1a1b1b046 100644 --- a/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp +++ b/src/QMCWaveFunctions/EinsplineSpinorSetBuilder.cpp @@ -55,8 +55,9 @@ std::unique_ptr EinsplineSpinorSetBuilder::createSPOSetFromXML(xmlNodePt { OhmmsAttributeSet a; + TinyVector TileFactor_do_not_use; a.add(H5FileName, "href"); - a.add(TileFactor, "tile"); + a.add(TileFactor_do_not_use, "tile", {}, TagStatus::DELETED); a.add(sortBands, "sort"); a.add(TileMatrix, "tilematrix"); a.add(twist_num_inp, "twistnum"); From d06fd1f81783d5f4aea038284ccc47ce3f7d2a4c Mon Sep 17 00:00:00 2001 From: "Paul R. C. Kent" Date: Fri, 25 Aug 2023 10:17:15 -0400 Subject: [PATCH 9/9] Set development version 3.17.9 --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index cc6225cbc8..271631d012 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,7 +15,7 @@ endif() ###################################################################### project( qmcpack - VERSION 3.17.1 + VERSION 3.17.9 LANGUAGES C CXX) # add the automatically determined parts of the RPATH