Merge pull request QMCPACK#5649 from PDoakORNL/eden_agreement_hamiltonian_changes

collect missing particle elements for energy density

Author: ye-luo

src/QMCHamiltonians/CoulombPBCAA.cpp

@@ -2,7 +2,7 @@
 // This file is distributed under the University of Illinois/NCSA Open Source License.
 // See LICENSE file in top directory for details.
 //
-// Copyright (c) 2022 QMCPACK developers.
+// Copyright (c) 2025 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
@@ -32,8 +32,6 @@ class OneDimCubicSplineLinearGrid;
 /** @ingroup hamiltonian
  *\brief Calculates the AA Coulomb potential using PBCs
  *
- * Functionally identical to CoulombPBCAA but uses a templated version of
- * LRHandler.
  */
 struct CoulombPBCAA : public OperatorBase, public ForceBase
 {
@@ -44,7 +42,7 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
   using OffloadSpline  = OneDimCubicSplineLinearGrid<LRCoulombSingleton::pRealType>;
 
   /// energy-optimized long range handle. Should be const LRHandlerType eventually
-  std::shared_ptr<LRHandlerType> AA;
+  std::shared_ptr<LRHandlerType> lr_aa_;
   /// energy-optimized short range pair potential
   std::shared_ptr<const RadFunctorType> rVs;
   /// the same as rVs but can be used inside OpenMP offload regions
@@ -85,10 +83,31 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
   Array<TraceReal, 1>* V_sample;
   Array<TraceReal, 1> V_const;
 #endif
-  ParticleSet& Ps;
 
+  /** Needs to be a mutable reference and not a copy since ParticleSet
+   *  copy misses important state which causes later access
+   *  violations.
+   *  This is a bad bad smell.
+   */
+  ParticleSet& Ps;
 
-  /** constructor */
+  /** constructor
+   *  Sadly there is significant conditional behavior here
+   *  the constructor does the operator calculation for the static
+   *  and only for active==true are most of the methods other than
+   *  nops.
+   *
+   *  Consistent side effects
+   *  * local copies of most of of pset refs basic properties
+   *    input psets will not be checked for consistency later!
+   *  *
+   *  Conditional side effects include
+   *  When ComputeForces == true
+   *  * turnOnPerParticleSK for ref
+   *  * updateSource(ref)
+   *  When is_activate == false
+   *
+   */
   CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces, bool use_offload);
 
   ~CoulombPBCAA() override;
@@ -97,6 +116,13 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
 
   void resetTargetParticleSet(ParticleSet& P) override;
 
+  /** evaluate just one walker
+   *
+   *  This is still called for batch AA IonIon
+   *  we later need to make sure for Ions that we assign from value_
+   *  to the v_samples in multi walker reference.
+   *  AA IonIon evaluation.
+   */
   Return_t evaluate(ParticleSet& P) override;
 
   void mw_evaluate(const RefVectorWithLeader<OperatorBase>& o_list,
@@ -124,6 +150,11 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
                          TrialWaveFunction& psi,
                          ParticleSet::ParticlePos& hf_terms,
                          ParticleSet::ParticlePos& pulay_terms) override;
+
+  /**
+   *  calls  eval(LR|SR){withForces} and updates eS and eL updates new_value_
+   *  new_value_ doesn't seem to ever be used anywhere
+   */
   void updateSource(ParticleSet& s) override;
 
   /** Do nothing */
@@ -209,7 +240,13 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
 
   /** constructor code factored out
    */
-  void initBreakup(ParticleSet& P);
+  void initBreakup(ParticleSet& ref_pset);
+
+  /** initialize everything for a static particle set
+   *  for static particle set hamiltonians no significant calculation
+   *  is done after this.
+   */
+  void inactiveInitialization(ParticleSet& ref);
 
   /** Compute the const part of the per particle coulomb self interaction potential.
    *  \param[out]  pp_consts   constant values for the particles self interaction

src/QMCHamiltonians/CoulombPBCAA.h

@@ -17,6 +17,8 @@ namespace qmcplusplus
 template<typename T>
 void combinePerParticleEnergies(const CrowdEnergyValues<T>& cev_in, std::vector<Vector<T>>& values_out)
 {
+  if (cev_in.begin() == cev_in.end())
+    return;
   const auto num_walkers   = cev_in.begin()->second.size();
   const auto num_particles = cev_in.begin()->second[0].size();
 
@@ -50,5 +52,5 @@ template void combinePerParticleEnergies<std::complex<double>>(const CrowdEnergy
 template void combinePerParticleEnergies<std::complex<float>>(const CrowdEnergyValues<std::complex<float>>& cev_in,
                                                               std::vector<Vector<std::complex<float>>>& values_out);
 
-  
+
 } // namespace qmcplusplus

src/QMCHamiltonians/Listener.cpp

@@ -414,8 +414,8 @@ void NonLocalECPotential::mw_evaluateImpl(const RefVectorWithLeader<OperatorBase
           auto& vi_samples = O_leader.mw_res_handle_.getResource().vi_samples;
           // CAUTION! This may not be so simple in the future
           int iw = j;
-          ve_samples(iw, batch_list[j].get().electron_id) += pairpots[j];
-          vi_samples(iw, batch_list[j].get().ion_id) += pairpots[j];
+          ve_samples(iw, batch_list[j].get().electron_id) += 0.5 * pairpots[j];
+          vi_samples(iw, batch_list[j].get().ion_id) += 0.5 * pairpots[j];
         }
 
 #ifdef DEBUG_NLPP_BATCHED
@@ -453,7 +453,7 @@ void NonLocalECPotential::mw_evaluateImpl(const RefVectorWithLeader<OperatorBase
         listener.report(iw, O_leader.getName(), ve_sample);
 
       for (const ListenerVector<Real>& listener : listeners->ion_values)
-        listener.report(iw, O_leader.getName(), vi_sample);
+        listener.report(iw, ion_listener_operator_name, vi_sample);
     }
     ve_samples = 0.0;
     vi_samples = 0.0;

src/QMCHamiltonians/NonLocalECPotential.cpp

@@ -13,20 +13,6 @@
 namespace qmcplusplus
 {
 
-// See src/QMCHamiltonians/tests/test_hamiltonian_factory for parsing tests
-static constexpr const char* const hamiltonian_xml = R"(
-<hamiltonian name="h0" type="generic" target="e">
-  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
-</hamiltonian>
-  )";
-
-static constexpr const char* const hamiltonian_eeei_xml = R"(
-<hamiltonian name="h0" type="generic" target="e">
-  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
-  <pairpot type="coulomb" name="ElecIon" source="ion" target="e"/>
-</hamiltonian>
-  )";
-
 HamiltonianPool MinimalHamiltonianPool::make_hamWithEE(Communicate* comm,
                                                        ParticleSetPool& particle_pool,
                                                        WaveFunctionPool& wavefunction_pool)
@@ -55,5 +41,33 @@ HamiltonianPool MinimalHamiltonianPool::makeHamWithEEEI(Communicate* comm,
   return hpool;
 }
 
+HamiltonianPool MinimalHamiltonianPool::makeHamWithEEEIII(Communicate* comm,
+                                                          ParticleSetPool& particle_pool,
+                                                          WaveFunctionPool& wavefunction_pool)
+{
+  HamiltonianPool hpool(particle_pool, wavefunction_pool, comm);
+  Libxml2Document doc;
+  doc.parseFromString(hamiltonian_eeeiii_xml);
+
+  xmlNodePtr root = doc.getRoot();
+  hpool.put(root);
+
+  return hpool;
+}
+
+HamiltonianPool MinimalHamiltonianPool::makeHamWithEEEIPS(Communicate* comm,
+                                                          ParticleSetPool& particle_pool,
+                                                          WaveFunctionPool& wavefunction_pool)
+{
+  HamiltonianPool hpool(particle_pool, wavefunction_pool, comm);
+  Libxml2Document doc;
+  doc.parseFromString(hamiltonian_eeeips_xml);
+
+  xmlNodePtr root = doc.getRoot();
+  hpool.put(root);
+
+  return hpool;
+}
+
 
 } // namespace qmcplusplus

src/QMCHamiltonians/tests/MinimalHamiltonianPool.cpp

@@ -2,7 +2,7 @@
 // This file is distributed under the University of Illinois/NCSA Open Source License.
 // See LICENSE file in top directory for details.
 //
-// Copyright (c) 2019 QMCPACK developers.
+// Copyright (c) 2025 QMCPACK developers.
 //
 // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Laboratory
 //
@@ -21,13 +21,62 @@ namespace qmcplusplus
 {
 class MinimalHamiltonianPool
 {
+  // See src/QMCHamiltonians/tests/test_hamiltonian_factory for parsing tests
+  static constexpr const char* const hamiltonian_xml = R"(
+<hamiltonian name="h0" type="generic" target="e">
+  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
+</hamiltonian>
+  )";
+
+  static constexpr const char* const hamiltonian_eeei_xml = R"(
+<hamiltonian name="h0" type="generic" target="e">
+  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
+  <pairpot type="coulomb" name="ElecIon" source="ion" target="e"/>
+</hamiltonian>
+  )";
+
+  static constexpr const char* const hamiltonian_eeeiii_xml = R"(
+<hamiltonian name="h0" type="generic" target="e">
+  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
+  <pairpot type="coulomb" name="ElecIon" source="ion" target="e"/>
+  <pairpot type="coulomb" name="IonIon" source="ion" target="ion"/>
+</hamiltonian>
+  )";
+
+  static constexpr const char* const hamiltonian_eeeips_xml = R"(
+<hamiltonian name="h0" type="generic" target="e">
+  <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
+  <pairpot type="coulomb" name="ElecIon" source="ion" target="e"/>
+  <pairpot type="pseudo" name="PseudoPot" source="ion0" wavefunction="psi0" format="xml">
+     <pseudo elementType="C" href="C.BFD.xml"/>
+  </pairpot>
+</hamiltonian>
+  )";
+
 public:
+  /// make a HamitonianPool with a primary hamiltonian with Coulombic electron electron interaction
   static HamiltonianPool make_hamWithEE(Communicate* comm,
                                         ParticleSetPool& particle_pool,
                                         WaveFunctionPool& wavefunction_pool);
+  /** make a HamitonianPool with a primary hamiltonian with Coulombic electron
+   *  electron interaction and electron ion interaction
+   */
   static HamiltonianPool makeHamWithEEEI(Communicate* comm,
                                          ParticleSetPool& particle_pool,
                                          WaveFunctionPool& wavefunction_pool);
+  /** make a HamitonianPool with a primary hamiltonian with Coulombic electron
+  *   electron interaction and electron ion interaction and ion ion interaction.
+  */
+  static HamiltonianPool makeHamWithEEEIII(Communicate* comm,
+                                           ParticleSetPool& particle_pool,
+                                           WaveFunctionPool& wavefunction_pool);
+
+  /** make a HamitonianPool with a primary hamiltonian with Coulombic electron
+  *   electron interaction and electron ion interaction and a pseudo potential component
+  */
+  static HamiltonianPool makeHamWithEEEIPS(Communicate* comm,
+                                           ParticleSetPool& particle_pool,
+                                           WaveFunctionPool& wavefunction_pool);
 };
 
 } // namespace qmcplusplus

src/QMCHamiltonians/tests/MinimalHamiltonianPool.h

@@ -190,10 +190,12 @@ TEST_CASE("NonLocalECPotential", "[hamiltonian]")
   testing::TestNonLocalECPotential::evaluateImpl(nl_ecp, elec, false, true);
   const auto value = nl_ecp.getValue();
 
-  CHECK(std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0) == Approx(value));
-  CHECK(std::accumulate(local_pots2.begin(), local_pots2.begin() + local_pots2.cols(), 0.0) == Approx(value));
-  CHECK(std::accumulate(ion_pots.begin(), ion_pots.begin() + ion_pots.cols(), 0.0) == Approx(-2.1047365665));
-  CHECK(std::accumulate(ion_pots2.begin(), ion_pots2.begin() + ion_pots2.cols(), 0.0) == Approx(-2.1047365665));
+  double total_localpots = std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0);
+  total_localpots += std::accumulate(ion_pots.begin(), ion_pots.begin() + ion_pots.cols(), 0.0);
+  CHECK(total_localpots == Approx(value));
+  double total_localpots2 = std::accumulate(local_pots[1], local_pots[1] + local_pots.cols(), 0.0);
+  total_localpots2 += std::accumulate(ion_pots[1], ion_pots[1] + ion_pots.cols(), 0.0);
+  CHECK(total_localpots2 == Approx(value));
 
   CHECK(!testing::TestNonLocalECPotential::didGridChange(nl_ecp));
 
@@ -203,16 +205,22 @@ TEST_CASE("NonLocalECPotential", "[hamiltonian]")
   testing::TestNonLocalECPotential::mw_evaluateImpl(nl_ecp, o_list, twf_list, p_list, false, listener_opt, true);
 
   CHECK(!testing::TestNonLocalECPotential::didGridChange(nl_ecp));
-  auto value2 = o_list[0].evaluateDeterministic(elec);
+  auto value_moved = o_list[0].evaluateDeterministic(elec);
 
-  CHECK(std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0) == Approx(value2));
+  total_localpots = std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0);
+  total_localpots += std::accumulate(ion_pots.begin(), ion_pots.begin() + ion_pots.cols(), 0.0);
+  CHECK(total_localpots == Approx(value_moved));
   // check the second walker which will be unchanged.
-  CHECK(std::accumulate(local_pots2[1], local_pots2[1] + local_pots2.cols(), 0.0) == Approx(value));
+  total_localpots2 = std::accumulate(local_pots[1], local_pots[1] + local_pots.cols(), 0.0);
+  total_localpots2 += std::accumulate(ion_pots[1], ion_pots[1] + ion_pots.cols(), 0.0);
+  CHECK(total_localpots2 == Approx(value));
 
-  // Randomizing grid does nothing for Na pp
   testing::TestNonLocalECPotential::mw_evaluateImpl(nl_ecp, o_list, twf_list, p_list, false, listener_opt, false);
-  auto value3 = o_list[0].evaluateDeterministic(p_list[0]);
-  CHECK(std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0) == Approx(value3));
+  auto value3     = o_list[0].evaluateDeterministic(p_list[0]);
+  total_localpots = std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0);
+  total_localpots += std::accumulate(ion_pots.begin(), ion_pots.begin() + ion_pots.cols(), 0.0);
+
+  CHECK(total_localpots == Approx(value3));
 }
 
 } // namespace qmcplusplus