[Codescan] Reject unsupported special fractional-derivative parameters instead of returning zero

[Stardust0831]

Global scan provenance

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

Problem

The fractional-derivative setup accepts any alpha <= p with p equal to 0 or 1. It later tells the user that the special algorithm only supports (p=1, alpha=1) and (p=0, alpha=0), but selecting special does not validate that the previously entered pair is supported.

Multiwfn/ext/2F2.f90

Lines 40 to 67 in d2b770a

	case("get")
	call get_alpha_level()
	return
	case("set_p")
	write (*, *) "Input fractional order of derivative (alpha); it should be <= 1"
	read (*, *) alpha_
	write (*, *) "Input integer order of internal derivative (p); it can be 0 or 1"
	read (*, *) p
	if (alpha_ > p) then
	write (*, *) "alpha can not be larger than p"
	write (*, *) ""
	cycle
	end if
	if (p /= 0 .and. p /= 1) then
	write (*, *) "p can be only 0 or 1"
	write (*, *) ""
	cycle
	end if
	case("set_auto")
	write (*, *) "Input fractional order of derivative (alpha); it should be <= 1"
	read (*, *) alpha_
	p = ceiling(alpha_)
	if (p < 0) p = 0
	write (*, '(A,I0)') "Automatically chosen integer order of internal derivative (p) is ", p
	case default
	write (*, *) "Try again..."
	cycle
	end select

Multiwfn/ext/2F2.f90

Lines 70 to 93 in d2b770a

	write (*, *) "Type `exact` for enabling evaluation of general analytic form of fractional derivative"
	write (*, *) "Type `approx` for interpolation of fractional derivative"
	write (*, *) "Type `special` for enabling evaluation of specific analytic form of fractional derivative"
	write (*, *) "For special, only several pairs are available:"
	write (*, *) " p alpha"
	write (*, *) " 1 1.0"
	write (*, *) " 0 0.0"
	do while (.true.)
	read (*, *) tempstr
	call to_lower(tempstr)
	select case(trim(tempstr))
	case("exact")
	algorithm_ = GENERAL_EXACT_EVALUATION
	case("approx")
	algorithm_ = INTERPOLATION_EVALUATION
	call prepare_interpolation()
	case("special")
	algorithm_ = SPECIAL_EXACT_EVALUATION
	case default
	write (*, *) "Try again..."
	cycle
	end select
	exit
	end do

The evaluation routine appears intended to raise an error for unsupported pairs, but initializes GTO_fractional_integral to 0 and then tests whether it differs from -1:

Multiwfn/ext/2F2.f90

Lines 116 to 165 in d2b770a

	real(8) function GTO_fractional_integral(n, x)
	use bspline_sub_module, only: db1val
	real(8), parameter :: one = 1e0_8, zero = 0e0_8
	integer, intent(in) :: n
	real(8), intent(in) :: x
	real(8) :: top, down
	real(8) :: gtop, gdown
	real(8) :: val2F2
	integer :: iout, inbvx
	if (x > zero) then
	write (*, *) "x is ", x
	error stop "x greater than zero!"
	end if
	val2F2 = -one
	GTO_fractional_integral = 0e0_8
	if (algorithm_ == GENERAL_EXACT_EVALUATION) then
	top = real(n + 1, kind=8)
	down = top + real(p, kind=8) - alpha_
	gtop = gamma(top)
	gdown = gamma(down)
	val2F2 = hyp2F2(top, top + one, down, down + one, x)
	GTO_fractional_integral = gtop / gdown * val2F2
	else if (algorithm_ == INTERPOLATION_EVALUATION) then
	top = real(n + 1, kind=8)
	down = top + real(p, kind=8) - alpha_
	gtop = gamma(top)
	gdown = gamma(down)
	inbvx = 1
	! normally, x shall be passed, but here some logic is broken since Bspline expects increasing values of x
	call db1val(-x, 0, t_2F2, n_2F2, k_2F2, Bs_2F2(:, n), val2F2, iout, inbvx)
	GTO_fractional_integral = gtop / gdown * val2F2
	else if (algorithm_ == SPECIAL_EXACT_EVALUATION) then
	select case(p)
	case(1)
	if (alpha_ == one) then
	GTO_fractional_integral = exp(x)
	end if
	case(0)
	if (alpha_ == zero) then
	GTO_fractional_integral = exp(x)
	end if
	end select
	if(GTO_fractional_integral /= -one) return
	write (*, *) "alpha: ", alpha_, " p: ", p
	error stop "No specialization for given alpha and p"
	else
	write (*, '(A,I0)') "Algorithm internal value is ", algorithm_
	error stop "Algorithm is not implemented!"
	end if
	end function GTO_fractional_integral

For an unsupported pair such as p=1, alpha=0.5, neither supported branch assigns a value. The result remains 0, 0 /= -1 is true, and the function returns before the diagnostic and error stop statements.

Impact

The program silently reports zero fractional integrals and therefore zero or corrupted fractional orbital derivatives for an explicitly selectable parameter combination. This is a scientific correctness failure: the calculation completes with plausible-looking numeric output instead of rejecting an unsupported method.

Suggested direction

• Validate the (p, alpha) pair immediately when special is selected and keep prompting or abort with a clear error if unsupported.
• Independently make the evaluator robust by using an explicit found flag or initializing the result to the sentinel actually tested.
• Avoid exact floating-point equality for future non-integer supported values; represent supported combinations explicitly or compare with a documented tolerance.
• Add tests for both supported pairs and at least one unsupported pair such as p=1, alpha=0.5, asserting that unsupported input cannot produce a silent zero result.