Multiwfn forum

Alexey

Member

Registered: 2024-06-28

Posts: 15

promolecule

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 sphatmraddens

i 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

Last edited by Alexey (2024-07-20 19:46:27)

sobereva

Tian Lu (Multiwfn developer)

From: Beijing

Registered: 2017-09-11

Posts: 2,011

Website

Re: promolecule

Please notice "ispheratm" in settings.ini. If you set it to 0, atomic wavefunction files will not be automatically sphericalized.