From eb744ccf36c780c48ab78f01ee83a022a85374de Mon Sep 17 00:00:00 2001
From: edoapra <edoardo.apra@gmail.com>
Date: Wed, 18 Jan 2023 16:34:22 -0800
Subject: [PATCH] use dftd3 method to address zero small denominators

---
 src/nwdft/xc/xc_vdw_util.F | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/nwdft/xc/xc_vdw_util.F b/src/nwdft/xc/xc_vdw_util.F
index 7c6443e03d..29f76bc1d4 100644
--- a/src/nwdft/xc/xc_vdw_util.F
+++ b/src/nwdft/xc/xc_vdw_util.F
@@ -42,6 +42,7 @@ C>
       double precision cnj !< [Input] \f$\mathrm{CN}^B(R)\f$
       double precision top,bottom,dist,c6_ref
       double precision cna,cnb
+      double precision c6_mem,dist_save
       integer i,j
       double precision eps
       parameter (eps=1d-90)
@@ -49,6 +50,8 @@ c
       c6cn=0.0d0
       top=0.0d0
       bottom=0.0d0
+      c6_mem=-1.d+99
+      dist_save=1.0d99
       do i=1,maxcn(iat)
          do j=1,maxcn(jat)
             c6_ref=c6ab(iat,jat,i,j,1)
@@ -56,6 +59,10 @@ c
                cna=c6ab(iat,jat,i,j,2)
                cnb=c6ab(iat,jat,i,j,3)
                dist=(cna-cni)**2+(cnb-cnj)**2
+               if (dist.lt.dist_save) then
+                  dist_save=dist
+                  c6_mem=c6_ref
+               endif
                top=top+dexp(k3*dist)*c6_ref
                bottom=bottom+dexp(k3*dist)
             endif
@@ -65,8 +72,7 @@ c
       if (bottom.gt.eps) then
          c6cn=top/bottom
       else
-         c6cn=0.0d0
-         call errquit(' c6cn returning zero ',0,0)
+         c6cn=c6_mem
       endif
       return
       end