Merge pull request QMCPACK#4775 from ye-luo/msd-nan-protection

Additional zero protection in MSD.
dcyang · Oct 16, 2023 · d0d456e · d0d456e
2 parents 6f7eefe + e18a785
commit d0d456e
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 12 deletions.
diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp
@@ -515,9 +515,8 @@ void MultiDiracDeterminant::evaluateDetsForPtclMove(const ParticleSet& P, int ia
     psiMinv_temp = psiMinv;
     for (size_t i = 0; i < NumPtcls; i++)
       psiV_temp[i] = psiV[*(it++)];
-    auto ratio_old_ref_det = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
-    curRatio               = ratio_old_ref_det;
-    InverseUpdateByColumn(psiMinv_temp, psiV_temp, workV1, workV2, WorkingIndex, ratio_old_ref_det);
+    curRatio = DetRatioByColumn(psiMinv_temp, psiV_temp, WorkingIndex);
+    InverseUpdateByColumn(psiMinv_temp, psiV_temp, workV1, workV2, WorkingIndex, curRatio);
     for (size_t i = 0; i < NumOrbitals; i++)
       TpsiM(i, WorkingIndex) = psiV[i];
   }

diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp
@@ -413,6 +413,8 @@ void MultiDiracDeterminant::acceptMove(ParticleSet& P, int iat, bool safe_to_del
   const int WorkingIndex = iat - FirstIndex;
   assert(WorkingIndex >= 0 && WorkingIndex < LastIndex - FirstIndex);
   assert(P.isSpinor() == is_spinor_);
+  if (curRatio == ValueType(0))
+    throw std::runtime_error("MultiDiracDeterminant::acceptMove curRatio is zero! Report a bug.\n");
   log_value_ref_det_ += convertValueToLog(curRatio);
   curRatio = ValueType(1);
   switch (UpdateMode)

diff --git a/src/QMCWaveFunctions/Fermion/MultiSlaterDetTableMethod.cpp b/src/QMCWaveFunctions/Fermion/MultiSlaterDetTableMethod.cpp
@@ -446,7 +446,9 @@ WaveFunctionComponent::PsiValue MultiSlaterDetTableMethod::ratioGrad(ParticleSet
     new_psi_ratio_to_new_ref_det_ = evalGrad_impl_no_precompute(P, iat, true, dummy);
 
   const int det_id = getDetID(iat);
-  curRatio         = Dets[det_id]->getRefDetRatio() * new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
+  curRatio         = Dets[det_id]->getRefDetRatio();
+  if (curRatio != PsiValue(0))
+    curRatio *= new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
   grad_iat += dummy;
   return curRatio;
 }
@@ -467,7 +469,9 @@ WaveFunctionComponent::PsiValue MultiSlaterDetTableMethod::ratioGradWithSpin(Par
     new_psi_ratio_to_new_ref_det_ = evalGradWithSpin_impl_no_precompute(P, iat, true, dummy, spindummy);
 
   const int det_id = getDetID(iat);
-  curRatio         = Dets[det_id]->getRefDetRatio() * new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
+  curRatio         = Dets[det_id]->getRefDetRatio();
+  if (curRatio != PsiValue(0))
+    curRatio *= new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
   grad_iat += dummy;
   spingrad_iat += spindummy;
   return curRatio;
@@ -501,7 +505,10 @@ void MultiSlaterDetTableMethod::mw_ratioGrad(const RefVectorWithLeader<WaveFunct
     auto& det                         = WFC_list.getCastedElement<MultiSlaterDetTableMethod>(iw);
     det.new_psi_ratio_to_new_ref_det_ = psi_list[iw];
     grad_new[iw] += dummy[iw];
-    ratios[iw] = det.curRatio = det.Dets[det_id]->getRefDetRatio() * psi_list[iw] / det.psi_ratio_to_ref_det_;
+    det.curRatio = det.Dets[det_id]->getRefDetRatio();
+    if (det.curRatio != PsiValue(0))
+      det.curRatio *= psi_list[iw] / det.psi_ratio_to_ref_det_;
+    ratios[iw] = det.curRatio;
   }
 }
 
@@ -548,7 +555,9 @@ WaveFunctionComponent::PsiValue MultiSlaterDetTableMethod::ratio(ParticleSet& P,
   Dets[det_id]->evaluateDetsForPtclMove(P, iat);
 
   new_psi_ratio_to_new_ref_det_ = computeRatio_NewMultiDet_to_NewRefDet(det_id);
-  curRatio = Dets[det_id]->getRefDetRatio() * new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
+  curRatio                      = Dets[det_id]->getRefDetRatio();
+  if (curRatio != PsiValue(0))
+    curRatio *= new_psi_ratio_to_new_ref_det_ / psi_ratio_to_ref_det_;
   return curRatio;
 }
 
@@ -605,9 +614,7 @@ void MultiSlaterDetTableMethod::mw_calcRatio(const RefVectorWithLeader<WaveFunct
       PsiValue psi_local(0);
       PRAGMA_OFFLOAD("omp parallel for reduction(+ : psi_local)")
       for (size_t i = 0; i < ndets; i++)
-      {
         psi_local += det_value_ptr_list_ptr[iw][i] * C_otherDs_ptr_list_ptr[iw][i];
-      }
       psi_list_ptr[iw] = psi_local;
     }
   }
@@ -616,7 +623,10 @@ void MultiSlaterDetTableMethod::mw_calcRatio(const RefVectorWithLeader<WaveFunct
   {
     auto& det                         = WFC_list.getCastedElement<MultiSlaterDetTableMethod>(iw);
     det.new_psi_ratio_to_new_ref_det_ = psi_list[iw];
-    ratios[iw] = det.curRatio = det.Dets[det_id]->getRefDetRatio() * psi_list[iw] / det.psi_ratio_to_ref_det_;
+    det.curRatio                      = det.Dets[det_id]->getRefDetRatio();
+    if (det.curRatio != PsiValue(0))
+      det.curRatio *= psi_list[iw] / det.psi_ratio_to_ref_det_;
+    ratios[iw] = det.curRatio;
   }
 }
 
@@ -632,7 +642,9 @@ void MultiSlaterDetTableMethod::evaluateRatios(const VirtualParticleSet& VP, std
     const OffloadVector<ValueType>& detValues0 = Dets[det_id]->getNewRatiosToRefDet();
 
     PsiValue psiNew = computeRatio_NewMultiDet_to_NewRefDet(det_id);
-    ratios[iat]     = Dets[det_id]->getRefDetRatio() * psiNew / psi_ratio_to_ref_det_;
+    ratios[iat]     = Dets[det_id]->getRefDetRatio();
+    if (ratios[iat] != ValueType(0))
+      ratios[iat] *= psiNew / psi_ratio_to_ref_det_;
   }
 }
 
@@ -1068,7 +1080,9 @@ void MultiSlaterDetTableMethod::evaluateDerivRatios(const VirtualParticleSet& VP
 
     // calculate VP ratios
     PsiValue psiNew = computeRatio_NewMultiDet_to_NewRefDet(det_id);
-    ratios[iat]     = Dets[det_id]->getRefDetRatio() * psiNew / psi_ratio_to_ref_det_;
+    ratios[iat]     = Dets[det_id]->getRefDetRatio();
+    if (ratios[iat] != ValueType(0))
+      ratios[iat] *= psiNew / psi_ratio_to_ref_det_;
 
     // calculate VP ratios derivatives
     if (recalculate)