Multiwfn official website: //www.umsyar.com/multiwfn. Multiwfn forum in Chinese: http://bbs.keinsci.com/wfn
You are not logged in.
Pages:1
Hi everyone! Please recommend literature that reviews the topological analysis of the ELF, LOL and Laplacian functions. The meaning of the critical points (3;+1), (3;-1) and (3;+3) in the context of these functions is not entirely clear
Will a function for outputting a .wfn file with promolecular wavefunction be added someday so that only xyz and "atomwfn" can be used to generate this file, and not a .wfn file with real wavefunction, calculated in the same basis as the atoms in "atomwfn"?
Hello! i need to use aspherical atomic density to build promolecular density, but i dont understand how to do it, because function 100(hidden option) > 10 give spherical avaraged density even if i use atom.wfn file with aspherical density. i rewrote subroutine sphatmraddens to get values for all (r;theta;phi) points, but i dont understand how to use this data for building promolecular density. i attach my code and file density.txt for some atom which i get with the help of my code. should i create new lagintpol in (r,theta,phi) space to use my data in building promolecular density with calcprodens(x,y,z,0)? or what? help me please
subroutine sphatmraddens use defvar use util use functions implicit real*8 (a-h,o-z) ! Declare variables real*8, allocatable :: radx(:), rady(:), radz(:), radpos(:), density(:) integer :: ntheta, nphi real*8 :: theta, phi integer :: irad, itheta, iphi, idx real*8 :: rad, x, y, z integer :: ifinish, iprogstp, iprogcrit, nradpt integer :: itmp parameter (truncrho = 1D-8, rlow = 0D0, rhigh = 12D0, nradpt = 200, ntheta = 50, nphi = 100) ! Allocate arrays allocate(radx(nradpt), rady(nradpt), radz(nradpt), radpos(nradpt), density(nradpt*ntheta*nphi)) ifinish = 0 iprogstp = 20 iprogcrit = iprogstp write(*,*) "Calculating..." ! Parallel loop !$OMP PARALLEL DO SHARED(density, radpos, ifinish, iprogcrit) PRIVATE(irad, itheta, iphi, rad, x, y, z, theta, phi, idx) schedule(dynamic) NUM_THREADS(nthreads) do irad = 1, nradpt rad = rlow + (rhigh - rlow) * (irad - 1) / (nradpt - 1) radpos(irad) = rad do itheta = 0, ntheta-1 theta = pi * itheta / (ntheta - 1) ! theta ranges from 0 to pi do iphi = 0, nphi-1 phi = 2 * pi * iphi / nphi ! phi ranges from 0 to 2*pi x = rad * sin(theta) * cos(phi) y = rad * sin(theta) * sin(phi) z = rad * cos(theta) idx = (irad - 1) * ntheta * nphi + itheta * nphi + iphi + 1 ! Compute density density(idx) = fdens(x, y, z) end do end do ifinish = ifinish + 1 if (ifinish >= iprogcrit) then call showprog(ifinish, nradpt) iprogcrit = iprogcrit + iprogstp end if end do !$OMP END PARALLEL DO ! Output results open(10, file="density.txt", status="replace") itmp = 0 do irad = 1, nradpt do itheta = 0, ntheta-1 theta = pi * itheta / (ntheta - 1) do iphi = 0, nphi-1 phi = 2 * pi * iphi / nphi idx = (irad - 1) * ntheta * nphi + itheta * nphi + iphi + 1 if (density(idx) > truncrho) then itmp = itmp + 1 write(10, "(' r=', f12.6, ' theta=', f12.6, ' phi=', f12.6, ' density=', f25.10, 'D0')") radpos(irad), theta, phi, density(idx) end if end do end do end do close(10) write(*,*) "The result has been output to density.txt in the current folder" end subroutine sphatmraddensi get correct values for all (r,theta,phi) points with this code, but i need to construct density with this values (do not pay attention that the density values are small, in my case it should be so)
file:https://drive.google.com/file/d/136n0jN … sp=sharing
p.s. may be my way to build promolecular density with aspherical atomic densities is bad? is there other way to do it?
i need to use molecule.xyz input and then analyze promolecular density, so 1000>17 doesnt suit me because i need molecule.wfn file for it
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 ifthank 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 subroutineand 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 NaNI 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?
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?
а нельзя использовать функцию 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?
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)?
Pages:1