Skip to content

Commit 9d3834e

Browse files
ye-luoye-luo
authored
Merge pull request QMCPACK#5490 from PDoakORNL/sk_est_4
Batched Structure Factor Estimator and unit tests
2 parents 2ee0c75 + faa261b commit 9d3834e

9 files changed

Lines changed: 1228 additions & 11 deletions

‎src/Estimators/CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@ set(QMCEST_SRC
4242
NEReferencePoints.cpp
4343
NESpaceGrid.cpp
4444
EnergyDensityEstimator.cpp
45-
StructureFactorInput.cpp)
45+
StructureFactorInput.cpp
46+
StructureFactorEstimator.cpp)
4647

4748
####################################
4849
# create libqmcestimators

‎src/Estimators/StructureFactorEstimator.cpp‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
//////////////////////////////////////////////////////////////////////////////////////
2+
// This file is distributed under the University of Illinois/NCSA Open Source License.
3+
// See LICENSE file in top directory for details.
4+
//
5+
// Copyright (c) 2024 QMCPACK developers.
6+
//
7+
// File developed by: Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory
8+
//
9+
// Some code refactored from: QMCHamiltonian/{SkEstimator.cpp, SkAllEstimator.cpp}
10+
//////////////////////////////////////////////////////////////////////////////////////
11+
12+
#include "StructureFactorEstimator.h"
13+
#include "StructureFactorInput.h"
14+
#include "ParticleSet.h"
15+
#include <LongRange/StructFact.h>
16+
17+
namespace qmcplusplus
18+
{
19+
20+
StructureFactorEstimator::StructureFactorEstimator(const StructureFactorInput& sfi,
21+
const ParticleSet& pset_ions,
22+
const ParticleSet& pset_elec,
23+
DataLocality data_locality)
24+
: OperatorEstBase(data_locality),
25+
input_(sfi),
26+
elns_(pset_elec),
27+
elec_num_species_(elns_.getSpeciesSet().getTotalNum()),
28+
ions_(pset_ions),
29+
ion_num_species_(ions_.getSpeciesSet().getTotalNum())
30+
{
31+
my_name_ = "StructureFactorEstimator";
32+
33+
num_kpoints_ = pset_ions.getSimulationCell().getKLists().getNumK();
34+
kshell_offsets_ = pset_ions.getSimulationCell().getKLists().getKShell();
35+
int max_kshell = kshell_offsets_.size() - 1;
36+
37+
rhok_tot_r_.resize(num_kpoints_);
38+
rhok_tot_i_.resize(num_kpoints_);
39+
sfk_e_e_.resize(num_kpoints_);
40+
rhok_e_.resize(num_kpoints_);
41+
// Legacy comment, but I don't trust it.
42+
//for values, we are including e-e structure factor, and e-Ion. So a total of NumIonSpecies+1 structure factors.
43+
//+2 for the real and imaginary parts of rho_k^e
44+
//
45+
// skAll seems to be written for e-e sf + complex rho_k^e
46+
data_.resize(3 * num_kpoints_);
47+
kmags_.resize(max_kshell);
48+
one_over_degeneracy_kshell_.resize(max_kshell);
49+
for (int ks = 0; ks < max_kshell; ks++)
50+
{
51+
kmags_[ks] = std::sqrt(pset_elec.getSimulationCell().getKLists().getKSQWorking()[kshell_offsets_[ks]]);
52+
one_over_degeneracy_kshell_[ks] = 1.0 / static_cast<Real>(kshell_offsets_[ks + 1] - kshell_offsets_[ks]);
53+
};
54+
}
55+
56+
StructureFactorEstimator::StructureFactorEstimator(const StructureFactorEstimator& sfe, DataLocality dl)
57+
: qmcplusplus::StructureFactorEstimator(sfe)
58+
{
59+
data_locality_ = dl;
60+
}
61+
62+
void StructureFactorEstimator::accumulate(const RefVector<MCPWalker>& walkers,
63+
const RefVector<ParticleSet>& psets,
64+
const RefVector<TrialWaveFunction>& wfns,
65+
const RefVector<QMCHamiltonian>& hams,
66+
RandomBase<FullPrecReal>& rng)
67+
{
68+
for (int iw = 0; iw < walkers.size(); ++iw)
69+
{
70+
Real weight = walkers[iw].get().Weight;
71+
72+
//sum over species
73+
std::copy(psets[iw].get().getSK().rhok_r[0], psets[iw].get().getSK().rhok_r[0] + num_kpoints_, rhok_tot_r_.begin());
74+
std::copy(psets[iw].get().getSK().rhok_i[0], psets[iw].get().getSK().rhok_i[0] + num_kpoints_, rhok_tot_i_.begin());
75+
for (int i = 1; i < elec_num_species_; ++i)
76+
accumulate_elements(psets[iw].get().getSK().rhok_r[i], psets[iw].get().getSK().rhok_r[i] + num_kpoints_,
77+
rhok_tot_r_.begin());
78+
for (int i = 1; i < elec_num_species_; ++i)
79+
accumulate_elements(psets[iw].get().getSK().rhok_i[i], psets[iw].get().getSK().rhok_i[i] + num_kpoints_,
80+
rhok_tot_i_.begin());
81+
82+
for (int k = 0; k < num_kpoints_; k++)
83+
{
84+
sfk_e_e_[k] += weight * (rhok_tot_r_[k] * rhok_tot_r_[k] + rhok_tot_i_[k] * rhok_tot_i_[k]);
85+
rhok_e_[k] += weight * std::complex<Real>{rhok_tot_r_[k], rhok_tot_i_[k]};
86+
}
87+
}
88+
}
89+
90+
void StructureFactorEstimator::registerOperatorEstimator(hdf_archive& file)
91+
{
92+
hdf_path hdf_name{my_name_};
93+
hdf_path path_variables = hdf_name / std::string_view("kpoints");
94+
file.push(path_variables, true);
95+
// hdf_archive wants non const references, if that code was better
96+
// this would be unecessary
97+
file.write(const_cast<std::vector<KPt>&>(ions_.getSimulationCell().getKLists().getKptsCartWorking()), "value");
98+
file.pop();
99+
}
100+
101+
void StructureFactorEstimator::write(hdf_archive& file)
102+
{
103+
hdf_path hdf_name{my_name_};
104+
file.push(hdf_name);
105+
// this is call rhok_e_e in the output of the legacy, but that is just wrong it is |rhok_e_e_|^2
106+
file.write(sfk_e_e_, "sfk_e_e");
107+
file.write(rhok_e_, "rhok_e_");
108+
}
109+
110+
void StructureFactorEstimator::collect(const RefVector<OperatorEstBase>& type_erased_operator_estimators)
111+
{
112+
int num_crowds = type_erased_operator_estimators.size();
113+
for (OperatorEstBase& crowd_oeb : type_erased_operator_estimators)
114+
{
115+
StructureFactorEstimator& crowd_sfe = dynamic_cast<StructureFactorEstimator&>(crowd_oeb);
116+
this->sfk_e_e_ += crowd_sfe.sfk_e_e_;
117+
this->rhok_e_ += crowd_sfe.rhok_e_;
118+
}
119+
OperatorEstBase::collect(type_erased_operator_estimators);
120+
}
121+
122+
void StructureFactorEstimator::startBlock(int steps) {}
123+
124+
UPtr<OperatorEstBase> StructureFactorEstimator::spawnCrowdClone() const
125+
{
126+
UPtr<StructureFactorEstimator> spawn(std::make_unique<StructureFactorEstimator>(*this, data_locality_));
127+
return spawn;
128+
}
129+
130+
} // namespace qmcplusplus

‎src/Estimators/StructureFactorEstimator.h‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
//////////////////////////////////////////////////////////////////////////////////////
2+
// This file is distributed under the University of Illinois/NCSA Open Source License.
3+
// See LICENSE file in top directory for details.
4+
//
5+
// Copyright (c) 2024 QMCPACK developers.
6+
//
7+
// File developed by: Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory
8+
//
9+
// Some code refactored from: QMCHamiltonian/{SkEstimator.h, SkAllEstimator.h}
10+
//////////////////////////////////////////////////////////////////////////////////////
11+
12+
#ifndef QMCPLUSPLUS_STRUCTUREFACTORESTIMATOR_H
13+
#define QMCPLUSPLUS_STRUCTUREFACTORESTIMATOR_H
14+
15+
#include "OperatorEstBase.h"
16+
17+
namespace qmcplusplus
18+
{
19+
20+
namespace testing
21+
{
22+
class StructureFactorAccess;
23+
}
24+
25+
class StructureFactorInput;
26+
27+
class StructureFactorEstimator : public OperatorEstBase
28+
{
29+
public:
30+
using QMCT = QMCTraits;
31+
using Real = QMCT::RealType;
32+
using FullPrecReal = QMCT::FullPrecRealType;
33+
using KPt = TinyVector<Real, QMCTraits::DIM>;
34+
35+
StructureFactorEstimator(const StructureFactorInput& sfi,
36+
const ParticleSet& pset_ions,
37+
const ParticleSet& pset_elec,
38+
DataLocality data_locality = DataLocality::crowd);
39+
40+
StructureFactorEstimator(const StructureFactorEstimator& sfe, DataLocality dl);
41+
42+
/** accumulate 1 or more walkers of EnergyDensity samples
43+
*/
44+
void accumulate(const RefVector<MCPWalker>& walkers,
45+
const RefVector<ParticleSet>& psets,
46+
const RefVector<TrialWaveFunction>& wfns,
47+
const RefVector<QMCHamiltonian>& hams,
48+
RandomBase<FullPrecReal>& rng) override;
49+
50+
/** start block entry point
51+
*/
52+
void startBlock(int steps) override;
53+
54+
UPtr<OperatorEstBase> spawnCrowdClone() const override;
55+
56+
void registerOperatorEstimator(hdf_archive& file) override;
57+
void write(hdf_archive& file) override;
58+
void collect(const RefVector<OperatorEstBase>& type_erased_operator_estimators) override;
59+
60+
long long getNumKPoints() const { return num_kpoints_; }
61+
const auto& getKLists() const { return ions_.getSimulationCell().getKLists(); };
62+
protected:
63+
// Testing functions
64+
const Vector<Real>& getSKElecElec() const { return sfk_e_e_; }
65+
const Vector<std::complex<Real>>& getRhoKElec() const { return rhok_e_; }
66+
67+
private:
68+
StructureFactorEstimator(const StructureFactorEstimator& obdm) = default;
69+
70+
const StructureFactorInput& input_;
71+
const ParticleSet& elns_;
72+
const int elec_num_species_;
73+
const ParticleSet& ions_;
74+
const int ion_num_species_;
75+
76+
/// number of k points
77+
long long num_kpoints_;
78+
79+
// std::vector<int> kshell_degeneracy_;
80+
/// kpts which belong to the ith-shell [kshell[i], kshell[i+1])
81+
std::vector<int> kshell_offsets_;
82+
83+
/// All the following are indexed by kshell
84+
/// instantaneous structure factor for a k shell
85+
std::vector<Real> kmags_;
86+
87+
std::vector<Real> one_over_degeneracy_kshell_;
88+
89+
/// Accumulated, its clearer to do it this way that use the base class data_ but means we don't use that base class infrastructure
90+
Vector<Real> sfk_e_e_;
91+
Vector<std::complex<Real>> rhok_e_;
92+
93+
/*@{
94+
* work area to reduce over electron species, they should probably just be complex as well.
95+
* but particle set stores the real and imaginary parts as two real vectors.
96+
*/
97+
Vector<Real> rhok_tot_r_;
98+
Vector<Real> rhok_tot_i_;
99+
/*@}*/
100+
101+
public:
102+
friend class testing::StructureFactorAccess;
103+
};
104+
105+
} // namespace qmcplusplus
106+
107+
#endif

‎src/Estimators/tests/CMakeLists.txt‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,11 +70,11 @@ set(SRCS
7070
test_EnergyDensityInput.cpp
7171
test_EnergyDensityEstimator.cpp
7272
test_StructureFactorInput.cpp
73+
test_StructureFactorEstimator.cpp
7374
)
7475

7576
add_executable(${UTEST_EXE} ${SRCS})
76-
target_link_libraries(${UTEST_EXE} catch_main qmcutil qmcestimators
77-
test_estimators_help utilities_for_test sposets_for_testing minimal_testing_pools)
77+
target_link_libraries(${UTEST_EXE} catch_main qmcutil qmcestimators test_estimators_help utilities_for_test sposets_for_testing minimal_testing_pools)
7878
if(USE_OBJECT_TARGET)
7979
target_link_libraries(
8080
${UTEST_EXE}

0 commit comments

Comments
 (0)