Fortran存储非球形势多极系数二进制文件运行TDSE报错咨询

背景 #

我正使用Fortran代码数值求解含时薛定谔方程（TDSE），运行代码需提供以无格式二进制格式存储的势多极系数，存储规则见对应示意图。
此前处理球形势时运行完全正常，我将所有数据存为单列，使用如下代码生成二进制文件：

program ar_pot
integer(kind=4) :: lmax, nr
real(kind=8) :: rmax
integer(kind=1) :: spherical, linear, even
integer :: i
real(kind=8) :: pi, r, V00
nr = 1024
lmax = 40
rmax = 160d0
spherical = 1
linear = 0
even = 0
pi = acos(-1d0)
open(1,file="potential",form="unformatted")
write(1) lmax,nr,rmax,spherical,linear,even !creating header
do i = 1, nr
    read(*,*) V00
    write(1) V00*sqrt(4.0*pi)
    write(*,*) i, V00
end do
close(1)
end program ar_pot

我通过./a.out <data.txt的方式读取含nr行单列数据的txt文件生成目标文件。

问题描述 #

目前我需要处理非球形势，需按照示意图中通用（otherwise）规则存储数据，但使用生成的非球形势、非线性、非偶对称势的二进制文件运行TDSE代码时，出现报错：forrtl: severe (151): allocatable array is already allocated，我需要排查是二进制文件存储逻辑错误还是程序本身故障。
为简化测试，我构造了一个仅第一列与原data.txt数据一致、其余列全为0的数组（列数为(lmax+1)^2，数据点以制表符分隔，注：header中存储的l_max与示意图中的l_max不一定一致），预期运行结果应与原球形势结果一致。我编写的生成对应二进制文件的代码如下：

program h_pot
integer(kind=4) :: lmax, nr
real(kind=8) :: rmax
integer(kind=1) :: spherical, linear, even
integer :: i
real(kind=8) :: pi, r
character*100000000 :: V00
nr = 1024
lmax = 40
rmax = 160d0
spherical = 0
linear = 0
even = 0
pi = acos(-1d0)
open(1,file="potential2",form="unformatted")
write(1) lmax,nr,rmax,spherical,linear,even
do i = 1, nr
    read(*,'(A)') V00
    write(1) trim(V00)
end do
close(1)
end program h_pot

我通过./a.out <data2的方式运行程序，data2为符合上述要求的txt文件，运行TDSE代码时返回的完整报错信息如下：

forrtl: severe (151): allocatable array is already allocated
Image PC Routine Line Source
libifcoremt.so.5 00007F081093280F for_alloc_allocat Unknown Unknown
tdse3d-splitop 000000000041F641 Unknown Unknown Unknown
tdse3d-splitop 000000000040CFB3 Unknown Unknown Unknown
tdse3d-splitop 0000000000402662 Unknown Unknown Unknown
libc-2.17.so 00007F080DC20B15 __libc_start_main Unknown Unknown
tdse3d-splitop 0000000000402559 Unknown Unknown Unknown

我对Fortran以及二进制格式数据存储的相关操作经验较少，请问上述代码存在哪些明显缺陷？

问题解答 #

现有代码的核心缺陷 #

• 写入数据类型完全错误：之前球形势版本写入二进制文件的是real(kind=8)类型的数值，而你新版代码直接将整行字符串裁剪后写入二进制，TDSE程序读取时预期解析双精度浮点数，实际读到的是字符串的二进制编码，会得到完全非法的数值，后续用错误数值作为数组长度申请内存时，就可能触发重复分配、长度非法等报错。
• 未正确解析多极系数：你需要将每行的(lmax+1)^2个浮点数从文本中逐个读出，转换为双精度浮点类型后再写入二进制文件，不能直接写入原始文本字符串。
• 不必要的大内存占用：声明长度为1亿的字符变量完全没有必要，还会占用过量内存，实际仅需足够存储一行数据的长度即可。

修正后的参考代码 #

program h_pot
    implicit none
    integer(kind=4) :: lmax, nr, n_coeff, i
    real(kind=8) :: rmax
    integer(kind=1) :: spherical, linear, even
    real(kind=8) :: pi
    real(kind=8), allocatable :: v_coeff(:)

    ! 参数配置和原有逻辑保持一致
    nr = 1024
    lmax = 40
    rmax = 160d0
    spherical = 0
    linear = 0
    even = 0
    pi = acos(-1d0)
    ! 计算单r点对应的多极系数总个数
    n_coeff = (lmax + 1) ** 2
    allocate(v_coeff(n_coeff))

    open(1, file="potential2", form="unformatted")
    ! 写入头信息，格式和原有逻辑完全对齐
    write(1) lmax, nr, rmax, spherical, linear, even

    do i = 1, nr
        ! 读取一行内的所有多极系数
        read(*, *) v_coeff
        ! 和旧版球形势逻辑对齐，第一列l=0,m=0的系数乘sqrt(4pi)做归一化
        v_coeff(1) = v_coeff(1) * sqrt(4.0d0 * pi)
        ! 写入该r点的所有双精度浮点系数
        write(1) v_coeff
    end do

    close(1)
    deallocate(v_coeff)
end program h_pot

后续排查建议 #

用修正后的代码重新生成二进制文件后再测试，如果仍出现分配报错，再排查TDSE主程序本身的内存分配逻辑问题，目前看二进制生成逻辑错误是首要诱因。

内容的提问来源于stack exchange，提问作者hh25