diff --git a/src/Particle/ParticleSetT.h b/src/Particle/ParticleSetT.h index 906e092adb..10b627696a 100644 --- a/src/Particle/ParticleSetT.h +++ b/src/Particle/ParticleSetT.h @@ -21,8 +21,6 @@ #ifndef QMCPLUSPLUS_PARTICLESETT_H #define QMCPLUSPLUS_PARTICLESETT_H -#include - #include "DTModes.h" #include "DynamicCoordinatesT.h" #include "MCCoordsT.hpp" @@ -38,6 +36,8 @@ #include "Walker.h" #include "type_traits/template_types.hpp" +#include + namespace qmcplusplus { /// forward declarations @@ -696,6 +696,7 @@ class ParticleSetT : public OhmmsElementBase { myTwist = t; } + inline const SingleParticlePos& getTwist() const { diff --git a/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.cpp b/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.cpp index ce4bb5e8aa..2469d3c1d2 100644 --- a/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.cpp +++ b/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.cpp @@ -17,6 +17,7 @@ #include "SplineR2RT.h" +#include "CPU/BLAS.hpp" #include "Concurrency/OpenMP.h" #include "QMCWaveFunctions/BsplineFactory/contraction_helper.hpp" #include "spline2/MultiBsplineEval.hpp" @@ -125,17 +126,27 @@ SplineR2RT::applyRotation( std::copy_n(spl_coefs, coefs_tot_size, coef_copy_->begin()); } - // Apply rotation the dumb way b/c I can't get BLAS::gemm to work... - for (auto i = 0; i < BasisSetSize; i++) { - for (auto j = 0; j < this->OrbitalSetSize; j++) { - const auto cur_elem = Nsplines * i + j; - auto newval{0.}; - for (auto k = 0; k < this->OrbitalSetSize; k++) { - const auto index = i * Nsplines + k; - newval += (*coef_copy_)[index] * rot_mat[k][j]; + if constexpr (std::is_same_v) { + // Here, ST should be equal to ValueType, which will be double for R2R. + // Using BLAS to make things faster + BLAS::gemm('N', 'N', this->OrbitalSetSize, BasisSetSize, + this->OrbitalSetSize, ST(1.0), rot_mat.data(), this->OrbitalSetSize, + coef_copy_->data(), Nsplines, ST(0.0), spl_coefs, Nsplines); + } + else { + // Here, ST is float but ValueType is double for R2R. Due to issues with + // type conversions, just doing naive matrix multiplication in this case + // to not lose precision on rot_mat + for (IndexType i = 0; i < BasisSetSize; i++) + for (IndexType j = 0; j < this->OrbitalSetSize; j++) { + const auto cur_elem = Nsplines * i + j; + FullPrecValueType newval{0.}; + for (IndexType k = 0; k < this->OrbitalSetSize; k++) { + const auto index = i * Nsplines + k; + newval += (*coef_copy_)[index] * rot_mat[k][j]; + } + spl_coefs[cur_elem] = newval; } - spl_coefs[cur_elem] = newval; - } } } diff --git a/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.h b/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.h index ece156ac1a..1e2a841e13 100644 --- a/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.h +++ b/src/QMCWaveFunctions/BsplineFactory/SplineR2RT.h @@ -40,8 +40,12 @@ class SplineR2RT : public BsplineSetT using SplineType = typename bspline_traits::SplineType; using BCType = typename bspline_traits::BCType; using DataType = ST; + using RealType = typename SPOSetT::RealType; + using IndexType = typename SPOSetT::IndexType; + using FullPrecValueType = double; using PointType = TinyVector; using SingleSplineType = UBspline_3d_d; + // types for evaluation results using TT = typename BsplineSetT::ValueType; using GGGVector = typename BsplineSetT::GGGVector; @@ -55,8 +59,6 @@ class SplineR2RT : public BsplineSetT using hContainer_type = VectorSoaContainer; using ghContainer_type = VectorSoaContainer; - using RealType = typename SPOSetT::RealType; - private: bool IsGamma; ///\f$GGt=G^t G \f$, transformation for tensor in LatticeUnit to diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 05c1fe018b..78cfb90d62 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -150,11 +150,17 @@ if(OHMMS_DIM MATCHES 3) endif(HAVE_EINSPLINE) # plane wave SPO - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWBasis.cpp PlaneWave/PWBasisT.cpp PlaneWave/PWParameterSet.cpp PlaneWave/PWOrbitalBuilder.cpp) + set(FERMION_SRCS ${FERMION_SRCS} + PlaneWave/PWBasis.cpp + PlaneWave/PWBasisT.cpp + PlaneWave/PWOrbitalSetT.cpp + PlaneWave/PWRealOrbitalSetT.cpp + PlaneWave/PWParameterSet.cpp + PlaneWave/PWOrbitalBuilder.cpp) if(QMC_COMPLEX) - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp PlaneWave/PWOrbitalSetT.cpp) + set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp) else() - set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp PlaneWave/PWRealOrbitalSetT.cpp) + set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp) endif(QMC_COMPLEX) if(NOT QMC_COMPLEX) diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp index f48ea6348a..46157f9b28 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp @@ -514,19 +514,19 @@ EinsplineSetBuilderT::AnalyzeTwists2( } TargetPtcl.setTwist(superFracs[twist_num_]); -#ifndef QMC_COMPLEX - // Check to see if supercell twist is okay to use with real wave - // functions - for (int dim = 0; dim < OHMMS_DIM; dim++) { - double t = 2.0 * superFracs[twist_num_][dim]; - if (std::abs(t - round(t)) > MatchingTol * 100) { - app_error() - << "Cannot use this super twist with real wavefunctions.\n" - << "Please recompile with QMC_COMPLEX=1.\n"; - APP_ABORT("EinsplineSetBuilder::AnalyzeTwists2"); + if constexpr (!IsComplex_t{}()) { + // Check to see if supercell twist is okay to use with real wave + // functions + for (int dim = 0; dim < OHMMS_DIM; dim++) { + double t = 2.0 * superFracs[twist_num_][dim]; + if (std::abs(t - round(t)) > MatchingTol * 100) { + app_error() + << "Cannot use this super twist with real wavefunctions.\n" + << "Please recompile with QMC_COMPLEX=1.\n"; + APP_ABORT("EinsplineSetBuilder::AnalyzeTwists2"); + } } } -#endif // Now check to see that each supercell twist has the right twists // to tile the primitive cell orbitals. const int numTwistsNeeded = std::abs(det(TileMatrix)); @@ -574,78 +574,80 @@ EinsplineSetBuilderT::AnalyzeTwists2( IncludeTwists.push_back(superSets[twist_num_][i]); // Now, find out which twists are distinct DistinctTwists.clear(); -#ifndef QMC_COMPLEX - std::vector copyTwists; - for (int i = 0; i < IncludeTwists.size(); i++) { - int ti = IncludeTwists[i]; - 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 = primcell_kpoints[tj]; - PosType sum = twist_i + twist_j; - PosType diff = twist_i - twist_j; - if (TwistPair(twist_i, twist_j)) - distinct = false; + if constexpr (!IsComplex_t{}()) { + std::vector copyTwists; + for (int i = 0; i < IncludeTwists.size(); i++) { + int ti = IncludeTwists[i]; + 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 = primcell_kpoints[tj]; + PosType sum = twist_i + twist_j; + PosType diff = twist_i - twist_j; + if (TwistPair(twist_i, twist_j)) + distinct = false; + } + if (distinct) + DistinctTwists.push_back(ti); + else + copyTwists.push_back(ti); } - if (distinct) - DistinctTwists.push_back(ti); - else - copyTwists.push_back(ti); - } - // Now determine which distinct twists require two copies - MakeTwoCopies.resize(DistinctTwists.size()); - for (int i = 0; i < DistinctTwists.size(); i++) { - MakeTwoCopies[i] = false; - int ti = DistinctTwists[i]; - PosType twist_i = primcell_kpoints[ti]; - for (int j = 0; j < copyTwists.size(); j++) { - int tj = copyTwists[j]; - PosType twist_j = primcell_kpoints[tj]; - if (TwistPair(twist_i, twist_j)) - MakeTwoCopies[i] = true; + // Now determine which distinct twists require two copies + MakeTwoCopies.resize(DistinctTwists.size()); + for (int i = 0; i < DistinctTwists.size(); i++) { + MakeTwoCopies[i] = false; + int ti = DistinctTwists[i]; + PosType twist_i = primcell_kpoints[ti]; + for (int j = 0; j < copyTwists.size(); j++) { + int tj = copyTwists[j]; + PosType twist_j = primcell_kpoints[tj]; + if (TwistPair(twist_i, twist_j)) + MakeTwoCopies[i] = true; + } + if (this->myComm->rank() == 0) { + std::array buf; + int length = std::snprintf(buf.data(), buf.size(), + "Using %d copies of twist angle [%6.3f, %6.3f, %6.3f]\n", + MakeTwoCopies[i] ? 2 : 1, twist_i[0], twist_i[1], + twist_i[2]); + if (length < 0) + throw std::runtime_error("Error generating string"); + app_log() << std::string_view(buf.data(), length); + app_log().flush(); + } } - if (this->myComm->rank() == 0) { - std::array buf; - int length = std::snprintf(buf.data(), buf.size(), - "Using %d copies of twist angle [%6.3f, %6.3f, %6.3f]\n", - MakeTwoCopies[i] ? 2 : 1, twist_i[0], twist_i[1], twist_i[2]); - if (length < 0) - throw std::runtime_error("Error generating string"); - app_log() << std::string_view(buf.data(), length); - app_log().flush(); + // Find out if we can make real orbitals + use_real_splines_ = true; + for (int i = 0; i < DistinctTwists.size(); i++) { + int ti = DistinctTwists[i]; + 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) + use_real_splines_ = false; } + if (use_real_splines_ && (DistinctTwists.size() > 1)) { + app_log() << "***** Use of real orbitals is possible, but not " + "currently implemented\n" + << " with more than one twist angle.\n"; + use_real_splines_ = false; + } + if (use_real_splines_) + app_log() << "Using real splines.\n"; + else + app_log() << "Using complex splines.\n"; } - // Find out if we can make real orbitals - use_real_splines_ = true; - for (int i = 0; i < DistinctTwists.size(); i++) { - int ti = DistinctTwists[i]; - 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) - use_real_splines_ = false; - } - if (use_real_splines_ && (DistinctTwists.size() > 1)) { - app_log() << "***** Use of real orbitals is possible, but not " - "currently implemented\n" - << " with more than one twist angle.\n"; + else { + DistinctTwists.resize(IncludeTwists.size()); + MakeTwoCopies.resize(IncludeTwists.size()); + for (int i = 0; i < IncludeTwists.size(); i++) { + DistinctTwists[i] = IncludeTwists[i]; + MakeTwoCopies[i] = false; + } use_real_splines_ = false; } - if (use_real_splines_) - app_log() << "Using real splines.\n"; - else - app_log() << "Using complex splines.\n"; -#else - DistinctTwists.resize(IncludeTwists.size()); - MakeTwoCopies.resize(IncludeTwists.size()); - for (int i = 0; i < IncludeTwists.size(); i++) { - DistinctTwists[i] = IncludeTwists[i]; - MakeTwoCopies[i] = false; - } - use_real_splines_ = false; -#endif } template diff --git a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSet.h b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSet.h index 225033214b..5add827a86 100644 --- a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSet.h +++ b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSet.h @@ -15,8 +15,8 @@ /** @file PWOrbitalSet.h * @brief Definition of member functions of Plane-wave basis set */ -#ifndef QMCPLUSPLUS_PLANEWAVE_ORBITALSETT_BLAS_H -#define QMCPLUSPLUS_PLANEWAVE_ORBITALSETT_BLAS_H +#ifndef QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H +#define QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H #include "QMCWaveFunctions/PlaneWave/PWBasis.h" #include "QMCWaveFunctions/SPOSet.h" diff --git a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h index d4e13de966..9103a16ee2 100644 --- a/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h +++ b/src/QMCWaveFunctions/PlaneWave/PWOrbitalSetT.h @@ -18,8 +18,8 @@ /** @file PWOrbitalSetT.h * @brief Definition of member functions of Plane-wave basis set */ -#ifndef QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H -#define QMCPLUSPLUS_PLANEWAVE_ORBITALSET_BLAS_H +#ifndef QMCPLUSPLUS_PLANEWAVE_ORBITALSETT_BLAS_H +#define QMCPLUSPLUS_PLANEWAVE_ORBITALSETT_BLAS_H #include "CPU/BLAS.hpp" #include "QMCWaveFunctions/PlaneWave/PWBasisT.h" diff --git a/src/QMCWaveFunctions/RotatedSPOsT.cpp b/src/QMCWaveFunctions/RotatedSPOsT.cpp index dabdc282a9..1aa8af8ada 100644 --- a/src/QMCWaveFunctions/RotatedSPOsT.cpp +++ b/src/QMCWaveFunctions/RotatedSPOsT.cpp @@ -307,7 +307,6 @@ template void RotatedSPOsT::buildOptVariables(const size_t nel) { -#if !defined(QMC_COMPLEX) /* Only rebuild optimized variables if more after-rotation orbitals are * needed Consider ROHF, there is only one set of SPO for both spin up and * down Nup > Ndown. nel_major_ will be set Nup. @@ -332,7 +331,6 @@ RotatedSPOsT::buildOptVariables(const size_t nel) buildOptVariables(created_m_act_rot_inds, created_full_rot_inds); } -#endif } template @@ -340,7 +338,6 @@ void RotatedSPOsT::buildOptVariables( const RotationIndices& rotations, const RotationIndices& full_rotations) { -#if !defined(QMC_COMPLEX) const size_t nmo = Phi->getOrbitalSetSize(); // create active rotations @@ -419,7 +416,6 @@ RotatedSPOsT::buildOptVariables( param[i] = this->myVars[i]; apply_rotation(param, false); } -#endif } template @@ -858,33 +854,32 @@ RotatedSPOsT::evaluateDerivatives(ParticleSetT& P, // possibly replace wit BLAS calls for (int i = 0; i < nel; i++) for (int j = 0; j < nmo; j++) - Bbar(i, j) = d2psiM_all(i, j) + - 2.0 * dot(myG_J[i], dpsiM_all(i, j)) + + Bbar(i, j) = d2psiM_all(i, j) + 2 * dot(myG_J[i], dpsiM_all(i, j)) + myL_J[i] * psiM_all(i, j); //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PART2 - const T* const A(psiM_all.data()); - const T* const Ainv(psiM_inv.data()); - const T* const B(Bbar.data()); - ValueMatrix T_mat; + const ValueType* const A(psiM_all.data()); + const ValueType* const Ainv(psiM_inv.data()); + const ValueType* const B(Bbar.data()); + ValueMatrix t; ValueMatrix Y1; ValueMatrix Y2; ValueMatrix Y3; ValueMatrix Y4; - T_mat.resize(nel, nmo); + t.resize(nel, nmo); Y1.resize(nel, nel); Y2.resize(nel, nmo); Y3.resize(nel, nmo); Y4.resize(nel, nmo); - BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), A, nmo, Ainv, nel, T(0.0), - T_mat.data(), nmo); - BLAS::gemm('N', 'N', nel, nel, nel, T(1.0), B, nmo, Ainv, nel, T(0.0), - Y1.data(), nel); - BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), T_mat.data(), nmo, Y1.data(), - nel, T(0.0), Y2.data(), nmo); - BLAS::gemm('N', 'N', nmo, nel, nel, T(1.0), B, nmo, Ainv, nel, T(0.0), - Y3.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nel, ValueType(1.0), A, nmo, Ainv, nel, + ValueType(0.0), t.data(), nmo); + BLAS::gemm('N', 'N', nel, nel, nel, ValueType(1.0), B, nmo, Ainv, nel, + ValueType(0.0), Y1.data(), nel); + BLAS::gemm('N', 'N', nmo, nel, nel, ValueType(1.0), t.data(), nmo, + Y1.data(), nel, ValueType(0.0), Y2.data(), nmo); + BLAS::gemm('N', 'N', nmo, nel, nel, ValueType(1.0), B, nmo, Ainv, nel, + ValueType(0.0), Y3.data(), nmo); // possibly replace with BLAS call Y4 = Y3 - Y2; @@ -894,8 +889,8 @@ RotatedSPOsT::evaluateDerivatives(ParticleSetT& P, if (kk >= 0) { const int p = m_act_rot_inds.at(i).first; const int q = m_act_rot_inds.at(i).second; - dlogpsi[kk] += T_mat(p, q); - dhpsioverpsi[kk] += T(-0.5) * Y4(p, q); + dlogpsi[kk] += t(p, q); + dhpsioverpsi[kk] += ValueType(-0.5) * Y4(p, q); } } } diff --git a/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp b/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp index 24a5087f79..e5c04d205f 100644 --- a/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp +++ b/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp @@ -633,7 +633,8 @@ TEMPLATE_TEST_CASE( xmlNodePtr sposet_builder = xmlFirstElementChild(root); xmlNodePtr sposet_ptr = xmlFirstElementChild(sposet_builder); - EinsplineSetBuilderT einSet(elec, ptcl.getPool(), c, sposet_builder); + EinsplineSetBuilderT einSet( + elec, ptcl.getPool(), c, sposet_builder); auto spo = einSet.createSPOSetFromXML(sposet_ptr); REQUIRE(spo);