[Codescan] Stop the legacy MBIS path for elements beyond argon

[Stardust0831]

Global scan provenance

This issue was found by the Codex global repository scan of non-unit-test files in Stardust0831/Multiwfn, pinned to commit d2b770a531eb7039153198d04520b5992f15c49d.

Problem

The legacy mbis_frj initialization supports shells only through argon. When it encounters an atom with nuclear charge above 18, it prints a warning and exits only the atom-initialization loop:

Multiwfn/ext/frj.f90

Lines 328 to 361 in d2b770a

	! estimate initial shell values for Zeff and population, only up to Ar at present....
	! Znuc = 0 indicate a ghost center, initial as Zeff=1 and with a tiny population
	mshell=0
	do iatm=1,ncenter
	znuc = a(iatm)%charge
	if (znuc.eq.0.0d0) then
	zeff(1,iatm) = 1.0d0
	shellpop(1,iatm)=1.0d-3
	else
	zeff(1,iatm) = 2.0d0*znuc
	endif
	if (znuc.gt.0.0d0 .and. znuc.le.2.0d0) then
	mshell(iatm) = 1
	shellpop(1,iatm)=znuc
	endif
	if (znuc.gt.2.0d0 .and. znuc.le.10.0d0) then
	mshell(iatm) = 2
	zeff(2,iatm) = 2.0d0
	shellpop(1,iatm)=2.0d0
	shellpop(2,iatm)=znuc-2.0d0
	endif
	if (znuc.gt.10.0d0 .and. znuc.le.18.0d0) then
	mshell(iatm) = 3
	zeff(2,iatm) = 2.0d0*sqrt(znuc)
	zeff(3,iatm) = 2.0d0
	shellpop(1,iatm)=2.0d0
	shellpop(2,iatm)=8.0d0
	shellpop(3,iatm)=znuc-10.0d0
	endif
	if (znuc.gt.18.0d0) then
	write(*,*)' Sorry, no MBIS for atoms beyond Ar yet....'
	exit
	endif
	enddo

The calculation then continues into the MBIS iterations with mshell=0 for the unsupported atom and for any atoms after the early loop exit:

Multiwfn/ext/frj.f90

Lines 371 to 418 in d2b770a

	write(*,*)
	write(*,*) "Performing MBIS iterations to refine atomic spaces..."
	lastcharge=0.0d0
	!Cycle each atom to calculate their charges
	do icyc=1,maxcyc
	if (ioutmedchg==1) write(*,*)
	if (icyc==1) then
	write(*,"(' Cycle',i5)") icyc
	else
	write(*,"(' Cycle',i5,' Maximum change:',f10.6)") icyc,varmax
	end if
	! calculate shell weights for each atomic grid point with the current shellpop and zeff
	! calculate the new shell population
	shelltmp = 0.0d0
	do iatm=1,ncenter
	gridatm%x=gridatmorg%x+a(iatm)%x !Move quadrature point to actual position in molecule
	gridatm%y=gridatmorg%y+a(iatm)%y
	gridatm%z=gridatmorg%z+a(iatm)%z
	do ipt=1+iradcut*sphpot,ntotpot
	! do ipt=1,ntotpot
	wtatm=0.0d0
	wtot = 0.0d0
	do jatm=1,ncenter
	dx = gridatm(ipt)%x - a(jatm)%x
	dy = gridatm(ipt)%y - a(jatm)%y
	dz = gridatm(ipt)%z - a(jatm)%z
	dis = dx*dx + dy*dy + dz*dz
	dis = sqrt(dis)
	do kshell=1,mshell(jatm)
	znuc = zeff(kshell,jatm)
	znorm = (znuc**3)/(8.0d0*pi)
	qshell = shellpop(kshell,jatm)
	! note that we do not need to multiply with the grid weight, as wtatm will only be used as relative values
	tmp = qshell*znorm*exp(-znuc*dis)
	if (tmp.lt.dencut) tmp = 0.0d0
	wtatm(kshell,jatm) = wtatm(kshell,jatm) + tmp
	wtot = wtot + wtatm(kshell,jatm)
	end do
	end do
	tmpden = tmpdens(iatm,ipt)
	if (wtot.gt.0.0d0 .and. tmpden.gt.eps) then
	do jatm=1,ncenter
	do kshell=1,mshell(jatm)
	shelltmp(kshell,jatm) = shelltmp(kshell,jatm) + wtatm(kshell,jatm)*tmpden/wtot
	end do
	end do
	endif

The public MBIS menu rejects only elements beyond radon, and option -3 calls this legacy implementation without an argon-limit check:

Multiwfn/population.f90

Lines 4700 to 4711 in d2b770a

	real*8 tmpdens(radpot*sphpot,ncenter) !tmpdens(ipt,iatm) corresponds to contribution of iatm to molecular density at grid ipt, and meantime multiplied by single-center integration weight at that point
	real*8 atmdis2min(ncenter)
	integer mshell(ncenter) !Actual number of shells of atoms
	integer :: maxcyc=500,ioutmedchg=0,ioutshell=0,ignorefar=1
	real*8 :: crit=0.0001D0,eps=1D-14,dencut=1D-10
	if (any(a%index>86)) then
	write(*,*) "Error: MBIS for elements beyond Rn is not supported"
	write(*,*) "Press ENTER button to exit"
	read(*,*)
	return
	end if

Multiwfn/population.f90

Lines 4719 to 4744 in d2b770a

	do while(.true.)
	write(*,*)
	call menutitle("MBIS",15,2)
	if (ignorefar==1) write(*,*) "-4 Toggle if reducing cost by ignoring atoms far from grid, current: Yes"
	if (ignorefar==0) write(*,*) "-4 Toggle if reducing cost by ignoring atoms far from grid, current: No"
	if (imode==0) write(*,*) "-3 Enter frj implementation of MBIS code"
	if (ioutshell==1) write(*,*) "-2 Toggle if outputting population and width of shells, current: Yes"
	if (ioutshell==0) write(*,*) "-2 Toggle if outputting population and width of shells, current: No"
	if (ioutmedchg==1) write(*,*) "-1 Toggle if outputting atomic charges during iterations, current: Yes"
	if (ioutmedchg==0) write(*,*) "-1 Toggle if outputting atomic charges during iterations, current: No"
	write(*,*) "0 Return"
	write(*,*) "1 Start calculation!"
	write(*,"(a,i4)") " 2 Set the maximum number of iterations, current:",maxcyc
	write(*,"(a,f10.6)") " 3 Set convergence criterion of atomic charges, current:",crit
	read(*,*) isel
	if (isel==0) then
	return
	else if (isel==-4) then
	if (ignorefar==1) then
	ignorefar=0
	else
	ignorefar=1
	end if
	else if (isel==-3) then
	call mbis_frj
	else if (isel==-2) then

Impact

For molecules containing K through Rn, the legacy path continues with missing pro-atom shell densities. It can converge to invalid charges and shell parameters while having told the user only that support is unavailable. If the unsupported atom is not last, later otherwise-supported atoms are also left uninitialized because the initialization loop exited early.

Suggested direction

• Validate every atom before starting mbis_frj and return with a clear error if any element exceeds Ar.
• Do not use exit as an error path from only the initialization loop.
• If broader element support is intended, initialize all supported shell models explicitly before enabling the calculation.
• Add a regression using at least one K/Fe-containing system and verify that the legacy path refuses to calculate rather than emitting charges.