In fact, it is well-known that electron density at nuclear position is ill-defined, because there is a cusp of electron density at this point (see https://en.wikipedia.org/wiki/Kato_theorem).
]]> if (iele>=1) then
if (r>atmrhocut(iele)) cycle
call genatmraddens(iele,rhoarr,npt) !Extract spherically averaged radial density of corresponding element at specific grids
if (idohess==0) then
call lagintpol(atmradpos(1:npt),rhoarr(1:npt),npt,r,term,der1r,der2r,2)
else if (idohess==1) then
call lagintpol(atmradpos(1:npt),rhoarr(1:npt),npt,r,term,der1r,der2r,3)
end if
elerho=elerho+term
der1rdr=der1r/r
derx=derx+der1rdr*rx
dery=dery+der1rdr*ry
derz=derz+der1rdr*rz
if (idohess==1) then
tmpval=(der2r-der1rdr)/r2
dxx=dxx+der1rdr+tmpval*rx2
dyy=dyy+der1rdr+tmpval*ry2
dzz=dzz+der1rdr+tmpval*rz2
dxy=dxy+tmpval*rx*ry
dyz=dyz+tmpval*ry*rz
dxz=dxz+tmpval*rx*rz
end if
end if