Re: modify function

dummy@example.com (sobereva) — Tue, 09 Jul 2024 05:28:45 +0000

Derivative at nuclear position cannot be obtained. In the case of good promolecular density, for each element, Multiwfn has an array recording radial density of its isolated state, and Lagrangian interpolation is used to obtain the value and derivatives at any radial distance. At nuclear position, the interpolation is not feasible, because it is the starting position of radial density.

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).

modify function

dummy@example.com (Alexey) — Mon, 08 Jul 2024 16:30:00 +0000

i need to modify this part of 'calchessmat_prodens' subroutine to get good promolecular density and calculate its derivatives for all elements. the following part of code doesnt calculate derivatives in the nuclear positions (but value of promolecular density is good calculated), however, in non-nuclear positions the derivative is calculated normally, what should I write in the following code to get derivatives in nuclear positionsn? (i ask this because if i use "output prop in point" for the iuserfunc==-2 function then i get derivatives at nuclear positions, but in the case below (elerho from my new calchessmat_prodens) i don't get derivatives of elerho in nuclear positions, only value)

                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