From e18a7854bd15199acd47087e46e20d680d05bb3f Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Fri, 13 Oct 2023 18:11:48 -0500 Subject: [PATCH] Additional zero/NaN protection in MSD. When the proposed configuration is on the nodal surface. Reference determinant curRatio goes to zero. This move should be rejected eventually by Monte Carlo. However, to make the code proceed, some NaN protection is needed. When Dets[det_id]->getRefDetRatio() is zero, new_psi_ratio_to_new_ref_det_ the sum of all the dets with respect to the ref one is NaN because it is the sum of all the 1/0. So keeping curRatio zero should be healthy. --- .../Fermion/MultiDiracDeterminant.2.cpp | 5 ++- .../Fermion/MultiDiracDeterminant.cpp | 2 ++ .../Fermion/MultiSlaterDetTableMethod.cpp | 32 +++++++++++++------ 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp index 9f8ee6c62f..dff4132ded 100644 --- 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 index 48a43a7d93..901feee457 100644 --- 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 index d268cbf82a..affa2ab6d9 100644 --- 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(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(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& 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)