Multiwfn forum

Alexey

Member

Registered: 2024-06-28

Posts: 15

function.f90

could you please explain to me what is the difference between calchessmat_dens_promol and calchessmat_prodens? which function i need to use to calculate promolecular density with built-in spherical atomic densities and then analyze it (calculate gradient, lapl on promolecular density)?

Last edited by Alexey (2024-07-06 23:48:24)

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: function.f90

You should use calchessmat_prodens. Note that currently for element index <=18, the proatomic density in RDG original paper is used in this subroutine, while for others, the proatomic density constructed by me is used. The quality of the former is much poorer than the latter, but faster to evaluate. If you need good promolecular density, I suggest modifying this subroutine to use the latter for all elements.

calchessmat_dens_promol is used to calculate density and derivatives based on promolecular wavefunction combined from isolated atomic wavefunction. In contrast, calchessmat_prodens only requires atomic coordinates and element information.

Alexey

Member

Registered: 2024-06-28

Posts: 15

Re: function.f90

а нельзя использовать функцию calcprodens(x,y,z,0) для генерации хорошей промолекулярной плотности а потом использовать ее для анализа? And could you please tell me, щn the basis of what function is topological analysis done iuserfunc==-2(calcprodens) How are derivatives and the Hessian matrix calculated?

Last edited by Alexey (2024-07-07 10:39:53)

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: function.f90

а нельзя использовать функцию calcprodens(x,y,z,0) для генерации хорошей промолекулярной плотности а потом использовать ее для анализа? And could you please tell me, щn the basis of what function is topological analysis done iuserfunc==-2(calcprodens) How are derivatives and the Hessian matrix calculated?

Please fully speak in English, otherwise I cannot exactly understand your question.

Alexey

Member

Registered: 2024-06-28

Posts: 15

Re: function.f90

sorry))
is it possible to use calcprodens(x,y,z,0) to generate "good" promolecular density and then analyze it? And could you please tell me, what function/subroutine is used for topological analysis of iuserfunc==-2(calcprodens) How are derivatives and the Hessian matrix calculated?

Last edited by Alexey (2024-07-07 11:07:55)

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: function.f90

Yes, calcprodens(x,y,z,0) is exactly what you need.

In "real*8 function userfunc(x,y,z)", you can find

case (-2) !Promolecular density
    userfunc=calcprodens(x,y,z,0)

So, if you set iuserfunc=-2 and perform topology analysis for user-defined function, then it is equivalent to perform topology analysis on promolecular density. In this case, the derivatives are evaluated fully numerically.

The topology analysis module calls "subroutine gencalchessmat" to obtain needed derivatives (gradient and possibly Hessian). If you find there is no specific code in this subroutine for evaluating analytic derivatives, that means the derivatives will be evaluated numerically automatically.

Alexey

Member

Registered: 2024-06-28

Posts: 15

Re: function.f90

I create my own function that creates promolecular density (input file xyz) with calcprodens(x,y,z,0) and then analyze it. should i use subroutine gencalchessmat to calc promolgrad and hess?

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: function.f90

If your function doesn't have capability of calculating analytic gradient and Hessian, you do not need to modify gencalchessmat. The gradient and Hessian will be calculated numerically and automatically in gencalchessmat.

Alexey

Member

Registered: 2024-06-28

Posts: 15

Re: function.f90

thank you soooo much
but so, im following to your advice to modife calchessmat_prodens code to generate good promoldens
i did this:

(i changed 'if ele<=18' to 'if ele >=118)

subroutine calchessmat_prodens(xin,yin,zin,elerho,elegrad,elehess)
use util
real*8 elerho,xin,yin,zin
real*8,optional :: elegrad(3),elehess(3,3)
real*8 posarr(200),rhoarr(200),tvec(3)
elerho=0D0
derx=0D0
dery=0D0
derz=0D0
dxx=0D0
dyy=0D0
dzz=0D0
dxy=0D0
dyz=0D0
dxz=0D0
idohess=0
if (present(elehess)) idohess=1

call getpointcell(xin,yin,zin,ic,jc,kc)
do icell=ic-PBCnx,ic+PBCnx
? ? do jcell=jc-PBCny,jc+PBCny
? ? ? ? do kcell=kc-PBCnz,kc+PBCnz
? ? ? ? ? ? call tvec_PBC(icell,jcell,kcell,tvec)
? ? ? ? ? ? do i=1,nfragatm
? ?? ? ? ? ? ? ?iatm=fragatm(i)
? ?? ? ? ? ? ? ?iele=a(iatm)%index
? ?? ? ? ? ? ? ?!rx=a(iatm)%x+tvec(1)-xin !Wrong code, older than 2022-Sep-18
? ?? ? ? ? ? ? ?!ry=a(iatm)%y+tvec(2)-yin
? ?? ? ? ? ? ? ?!rz=a(iatm)%z+tvec(3)-zin
? ?? ? ? ? ? ? ?rx=xin-tvec(1)-a(iatm)%x !Relative x
? ?? ? ? ? ? ? ?ry=yin-tvec(2)-a(iatm)%y
? ?? ? ? ? ? ? ?rz=zin-tvec(3)-a(iatm)%z
? ?? ? ? ? ? ? ?rx2=rx*rx
? ?? ? ? ? ? ? ?ry2=ry*ry
? ?? ? ? ? ? ? ?rz2=rz*rz
? ?? ? ? ? ? ? ?r2=rx2+ry2+rz2
? ?? ? ? ? ? ? ?r=dsqrt(r2)
? ?? ? ? ? ? ? ?if (iele>=118) then !H~Ar, use Weitao Yang's fitted parameters as original RDG paper
? ?? ? ? ? ? ? ? ? ?if (atomdenscut==1) then !Tight cutoff, for CHNO corresponding to cutoff at rho=0.00001
? ? ? ?? ? ? ? ? ? ? ? ?if (iele==1.and.r2>25D0) then !H, 6.63^2=43.9569. But this seems to be unnecessarily large, so I use 5^2=25
? ? ? ? ? ?? ? ? ? ? ? ? ? ?cycle
? ? ? ?? ? ? ? ? ? ? ? ?else if (iele==6.and.r2>58.6756D0) then !C, 7.66^2=58.6756
? ? ? ? ? ?? ? ? ? ? ? ? ? ?cycle
? ? ? ?? ? ? ? ? ? ? ? ?else if (iele==7.and.r2>43.917129D0) then !N, 6.627^2=43.917129
? ? ? ? ? ?? ? ? ? ? ? ? ? ?cycle
? ? ? ?? ? ? ? ? ? ? ? ?else if (iele==8.and.r2>34.9281D0) then !O, 5.91^2=34.9281
? ? ? ? ? ?? ? ? ? ? ? ? ? ?cycle
? ? ? ?? ? ? ? ? ? ? ? ?else if (r2>(2.5D0*vdwr(iele))**2) then !Other cases, larger than 2.5 times of its vdw radius will be skipped
? ? ? ? ? ?? ? ? ? ? ? ? ? ?cycle
? ? ? ?? ? ? ? ? ? ? ? ?end if
? ?? ? ? ? ? ? ? ? ?else if (atomdenscut==2) then !Medium cutoff, the result may be not as accurate as atomdenscut==1, but much more cheaper
? ? ? ?? ? ? ? ? ? ? ? ?if (r2>(2.2D0*vdwr(iele))**2) cycle
? ?? ? ? ? ? ? ? ? ?else if (atomdenscut==3) then !Loose cutoff, the most inaccurate
? ? ? ?? ? ? ? ? ? ? ? ?if (r2>(1.8D0*vdwr(iele))**2) cycle
? ?? ? ? ? ? ? ? ? ?else if (atomdenscut==4) then !Foolish cutoff, you need to know what you are doing
? ? ? ?? ? ? ? ? ? ? ? ?if (r2>(1.5D0*vdwr(iele))**2) cycle
? ?? ? ? ? ? ? ? ? ?end if
? ? ? ?? ? ? ? ? ? ?r2_1d5=r2**1.5D0
? ? ? ?? ? ? ? ? ? ?do iSTO=1,3
? ? ? ? ? ?? ? ? ? ? ? ?if (YWTatomcoeff(iele,iSTO)==0D0) cycle
? ? ? ? ? ?? ? ? ? ? ? ?expterm=YWTatomexp(iele,iSTO)
? ? ? ? ? ?? ? ? ? ? ? ?term=YWTatomcoeff(iele,iSTO)*dexp(-r/expterm)
? ? ? ? ? ?? ? ? ? ? ? ?elerho=elerho+term
? ? ? ? ? ?? ? ? ? ? ? ?if (r==0D0) cycle !Derivative of STO at nuclei is pointless
? ? ? ? ? ?? ? ? ? ? ? ?tmp=term/expterm/r
? ? ? ? ? ?? ? ? ? ? ? ?derx=derx-tmp*rx !Calculating gradient doesn't cost detectable time, so always calculate it
? ? ? ? ? ?? ? ? ? ? ? ?dery=dery-tmp*ry
? ? ? ? ? ?? ? ? ? ? ? ?derz=derz-tmp*rz
? ? ? ? ? ?? ? ? ? ? ? ?if (idohess==1) then
? ? ? ? ? ? ? ?? ? ? ? ? ? ?tmp1=1/r2_1d5/expterm
? ? ? ? ? ? ? ?? ? ? ? ? ? ?tmp2=1/r2/(expterm*expterm)
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dxx=dxx+term*(tmp1*rx2-1/r/expterm+tmp2*rx2)
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dyy=dyy+term*(tmp1*ry2-1/r/expterm+tmp2*ry2)
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dzz=dzz+term*(tmp1*rz2-1/r/expterm+tmp2*rz2)
? ? ? ? ? ? ? ?? ? ? ? ? ? ?tmp=term*(tmp1+tmp2)
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dxy=dxy+rx*ry*tmp
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dyz=dyz+ry*rz*tmp
? ? ? ? ? ? ? ?? ? ? ? ? ? ?dxz=dxz+rx*rz*tmp
? ? ? ? ? ?? ? ? ? ? ? ?end if
? ? ? ?? ? ? ? ? ? ?end do
? ?? ? ? ? ? ? ?else !Heavier than Ar
? ? ? ? ? ? ? ? ? ? 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 !See promolecular_grid routine in props.f90 of NCIplot
? ? ? ? ? ?? ? ? ? ? ? ?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
? ? ? ? ? ? end do
? ? ? ? end do
? ? end do
end do
if (present(elegrad)) then
? ? elegrad(1)=derx
? ? elegrad(2)=dery
? ? elegrad(3)=derz
end if
if (idohess==1) then
? ? elehess(1,1)=dxx
? ? elehess(2,2)=dyy
? ? elehess(3,3)=dzz
? ? elehess(1,2)=dxy
? ? elehess(2,3)=dyz
? ? elehess(1,3)=dxz
? ? elehess(2,1)=dxy
? ? elehess(3,2)=dyz
? ? elehess(3,1)=dxz
end if
end subroutine

and if i try to "Output all properties at a point" (point is the O atom in H2O) i get follow for promolecular density (promolecular density is good calculated, but its derivatives are not calculated)

Density of electrons:? 0.3441456709E+00
Reduced density gradient:? 0.1000000000E+03
Note: Matrix diagonalization exceed max cycle before convergence
Sign(lambda2)*rho:? ? ? ? ? ? ? ?NaN
ESP from nuclear charges:? 0.1000000000E+04
van der Waals potential (probe atom: C ):? 0.1280973043+126 kcal/mol
User-defined real space function:? ? ? ? ? ? ? ?NaN

Note: Below information are for electron density

Components of gradient in x/y/z are:
? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN
Norm of gradient is:? ? ? ? ? ? ? ?NaN

Components of Laplacian in x/y/z are:
? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN
Total:? ? ? ? ? ? ? ?NaN

Hessian matrix:
? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN
? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN
? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN? ? ? ? ? ? ? ?NaN

Last edited by Alexey (2024-07-07 22:45:48)

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: function.f90

You can see

Note: Below information are for electron density

If currently wavefunction information is not available, then the data must be NaN.

If you want to obtain derivative information for e.g. real space function 100 (current user-defined function), you can enter main function 1, input f100, then input the coordinate of interest, and then various data at this point will be printed, and derivative information (like gradient and Hessian) will also be printed for the real space function 100.