From 2760da35c40ef78f737a8b3f0890482bf699474a Mon Sep 17 00:00:00 2001 From: Mingyu Zhu Date: Sun, 31 May 2026 20:43:18 +0800 Subject: [PATCH 1/4] enhance efficiency of `rfac.read_xx` for large files --- python/rfac.py | 1139 +++++++++++++++++++++++++----------------------- 1 file changed, 590 insertions(+), 549 deletions(-) diff --git a/python/rfac.py b/python/rfac.py index 0e244bc0..e7d665a7 100644 --- a/python/rfac.py +++ b/python/rfac.py @@ -26,7 +26,7 @@ from pfac import fac from pfac import const from pfac import util -import os +import os, datetime from multiprocessing import Pool, cpu_count def e2v(e, m=0): @@ -400,24 +400,22 @@ def get_lname(line): return get_lcomplex, get_lsname, get_lname -def _read_value(lines, cls): +def _read_value(lines, idx, cls): """ - Pop one line from lines, - split line by '=' and cast the value, then rturn the value and the rest of - lines. + Split the line of `idx` by '=' and cast the value, + then return the value and the moved index (`idx+1`). """ - vals = lines[0].split('=') + vals = lines[idx].split('=') if len(vals) > 1: val = cls(vals[1].strip()) else: val = '' - return val, lines[1:] + return val, idx+1 - -def _get_header(lines): - """ Read fac header """ +def _get_header(lines, idx=0): + """Read fac file header.""" header = {} - a = lines[0][4:-1].split('[') + a = lines[idx][4:-1].split('[') header['FAC'] = a[0] if len(a) > 1: a = a[1][:-1].split('.') @@ -427,389 +425,438 @@ def _get_header(lines): else: header['nthreads'] = 0 header['uta'] = 0 - header['utaci'] = 0 - lines = lines[1:] - header['Endian'], lines = _read_value(lines, int) - header['TSess'], lines = _read_value(lines, int) - header['Type'], lines = _read_value(lines, int) - header['Verbose'], lines = _read_value(lines, int) - key = lines[0].split('\t')[0] - header[key], lines = _read_value(lines, float) + header['utaci'] = 0 + idx += 1 + header['Endian'], idx = _read_value(lines, idx, int) + header['TSess'], idx = _read_value(lines, idx, int) + header['Type'], idx = _read_value(lines, idx, int) + header['Verbose'], idx = _read_value(lines, idx, int) + key = lines[idx].split('\t')[0] + header[key], idx = _read_value(lines, idx, float) a = key.split(' ') header['asym'] = a[0] header['Z'] = int(header[key]) - header['NBlocks'], lines = _read_value(lines, int) - return header, lines - + header['NBlocks'], idx = _read_value(lines, idx, int) + return header, idx def read_lev(filename): """ read *a.lev / *a.en file. """ with open(filename, 'r') as f: lines = f.readlines() - # header - header, lines = _get_header(lines) + # 1. parse header + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - ind, e0 = lines[0].split('=')[-1].split(',') + return header, () + + # 2. parse E0 + ind, e0 = lines[idx].split('=')[-1].split(',') header['E0_index'] = int(ind) header['E0'] = float(e0) - - lines = lines[2:] - if len(lines) > 3: - line0 = lines[3] + idx += 2 + if idx + 3 < len(lines): + line0 = lines[idx + 3] else: line0 = None get_lcomplex, get_lsname, get_lname = _wrap_get_length(line0) - def read_blocks(lines): + # 3. parse blocks + blocks = [] + while idx < len(lines): + # skip empty lines + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - nlev, lines = _read_value(lines, int) - # read the values - block['ILEV'] = np.zeros(nlev, dtype=int) - block['IBASE'] = np.zeros(nlev, dtype=int) - block['ENERGY'] = np.zeros(nlev, dtype=float) - block['P'] = np.zeros(nlev, dtype=int) - block['VNL'] = np.zeros(nlev, dtype=int) - block['2J'] = np.zeros(nlev, dtype=int) - block['ncomplex'] = np.chararray(nlev, itemsize=32) - block['sname'] = np.chararray(nlev, itemsize=48) - block['name'] = np.chararray(nlev, itemsize=128) - lines = lines[1:] - - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks - a = line.split() - block['ILEV'][i] = int(a[0]) - block['IBASE'][i] = int(a[1]) - block['ENERGY'][i] = float(a[2]) - block['P'][i] = int(a[3]) - block['VNL'][i] = int(a[4]) - block['2J'][i] = int(a[5]) - block['ncomplex'][i] = a[6] - block['sname'][i] = a[7] - block['name'][i] = a[8] - - return (block, ) - - return header, read_blocks(lines) + block['NELE'], idx = _read_value(lines, idx, int) + nlev, idx = _read_value(lines, idx, int) + idx += 1 + + ilev, ibase, energy, p, vnl, twoj = [], [], [], [], [], [] + ncomplex, sname, name = [], [], [] + + while idx < len(lines): + line = lines[idx] + if not line.strip(): # if empty + idx += 1 + break + a = line.split() + ilev.append(int(a[0])) + ibase.append(int(a[1])) + energy.append(float(a[2])) + p.append(int(a[3])) + vnl.append(int(a[4])) + twoj.append(int(a[5])) + ncomplex.append(a[6]) + sname.append(a[7]) + name.append(a[8]) + idx += 1 + + block['ILEV'] = np.array(ilev, dtype=int) + block['IBASE'] = np.array(ibase, dtype=int) + block['ENERGY'] = np.array(energy, dtype=float) + block['P'] = np.array(p, dtype=int) + block['VNL'] = np.array(vnl, dtype=int) + block['2J'] = np.array(twoj, dtype=int) + block['ncomplex'] = np.array(ncomplex, dtype='U32') + block['sname'] = np.array(sname, dtype='U48') + block['name'] = np.array(name, dtype='U128') + blocks.append(block) + + return header, tuple(blocks) def read_en(filename): return read_lev(filename) def read_enf(filename): - """ read en file with B&E """ + """ read en file with B&E. """ with open(filename, 'r') as f: lines = f.readlines() - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - nlev, lines = _read_value(lines, int) - block['EFIELD'], lines = _read_value(lines, float) - block['BFIELD'], lines = _read_value(lines, float) - block['FANGLE'], lines = _read_value(lines, float) - lines = lines[1:] - block['ilev'] = np.zeros(nlev, dtype=int) - block['energy'] = np.zeros(nlev, dtype=float) - block['pbasis'] = np.zeros(nlev, dtype=int) - block['mbasis'] = np.zeros(nlev, dtype=int) - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks + block['NELE'], idx = _read_value(lines, idx, int) + nlev, idx = _read_value(lines, idx, int) + block['EFIELD'], idx = _read_value(lines, idx, float) + block['BFIELD'], idx = _read_value(lines, idx, float) + block['FANGLE'], idx = _read_value(lines, idx, float) + idx += 1 + + ilev, energy, pbasis, mbasis = [], [], [], [] + while idx < len(lines): + line = lines[idx] + if not line.strip(): + idx += 1 + break a = line.split() - block['ilev'][i] = int(a[0]) - block['energy'][i] = float(a[1]) - block['pbasis'][i] = int(a[2]) - block['mbasis'][i] = int(a[3]) - return (block, ) + ilev.append(int(a[0])) + energy.append(float(a[1])) + pbasis.append(int(a[2])) + mbasis.append(int(a[3])) + idx += 1 + + block['ilev'] = np.array(ilev, dtype=int) + block['energy'] = np.array(energy, dtype=float) + block['pbasis'] = np.array(pbasis, dtype=int) + block['mbasis'] = np.array(mbasis, dtype=int) + blocks.append(block) + return header, tuple(blocks) - return header, read_blocks(lines) def read_trf(filename): with open(filename, 'r') as f: lines = f.readlines() - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['MULTIP'], lines = _read_value(lines, int) - block['GAUGE'], lines = _read_value(lines, int) - block['MODE'], lines = _read_value(lines, int) - block['EFIELD'], lines = _read_value(lines, float) - block['BFIELD'], lines = _read_value(lines, float) - block['FANGLE'], lines = _read_value(lines, float) - - block['upper_index'] = np.zeros(ntrans, dtype=int) - block['lower_index'] = np.zeros(ntrans, dtype=int) - block['upper_pbasis'] = np.zeros(ntrans, dtype=int) - block['lower_pbasis'] = np.zeros(ntrans, dtype=int) - block['upper_mbasis'] = np.zeros(ntrans, dtype=int) - block['lower_mbasis'] = np.zeros(ntrans, dtype=int) - block['energy'] = np.zeros(ntrans, dtype=float) - block['rate'] = np.zeros(ntrans, dtype=float) - nm = 2*abs(block['MULTIP'])+1 - block['mrate'] = np.zeros((ntrans,nm), dtype=float) + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['MULTIP'], idx = _read_value(lines, idx, int) + block['GAUGE'], idx = _read_value(lines, idx, int) + block['MODE'], idx = _read_value(lines, idx, int) + block['EFIELD'], idx = _read_value(lines, idx, float) + block['BFIELD'], idx = _read_value(lines, idx, float) + block['FANGLE'], idx = _read_value(lines, idx, float) + + nm = 2 * abs(block['MULTIP']) + 1 + upper_idx, lower_idx = [], [] + upper_pb, lower_pb = [], [] + upper_mb, lower_mb = [], [] + energy, rate = [], [] + mrate_list = [] + j = 0 - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks - im = i%nm - block['mrate'][j,im] = float(line[66:80]) - if im != nm-1: - continue - block['upper_index'][j] = int(line[:6]) - block['upper_pbasis'][j] = int(line[6:13]) - block['upper_mbasis'][j] = int(line[13:17]) - block['lower_index'][j] = int(line[17:24]) - block['lower_pbasis'][j] = int(line[24:31]) - block['lower_mbasis'][j] = int(line[31:35]) - block['energy'][j] = float(line[38:52]) - block['rate'][j] = float(line[94:108]) - j += 1 - - return (block, ) + current_mrate_row = [0.0] * nm + im_count = 0 + + while idx < len(lines): + line = lines[idx] + if not line.strip(): + idx += 1 + break + im = im_count % nm + current_mrate_row[im] = float(line[66:80]) + if im == nm - 1: + upper_idx.append(int(line[:6])) + upper_pb.append(int(line[6:13])) + upper_mb.append(int(line[13:17])) + lower_idx.append(int(line[17:24])) + lower_pb.append(int(line[24:31])) + lower_mb.append(int(line[31:35])) + energy.append(float(line[38:52])) + rate.append(float(line[94:108])) + mrate_list.append(current_mrate_row.copy()) + current_mrate_row = [0.0] * nm + j += 1 + im_count += 1 + idx += 1 + block['upper_index'] = np.array(upper_idx, dtype=int) + block['lower_index'] = np.array(lower_idx, dtype=int) + block['upper_pbasis'] = np.array(upper_pb, dtype=int) + block['lower_pbasis'] = np.array(lower_pb, dtype=int) + block['upper_mbasis'] = np.array(upper_mb, dtype=int) + block['lower_mbasis'] = np.array(lower_mb, dtype=int) + block['energy'] = np.array(energy, dtype=float) + block['rate'] = np.array(rate, dtype=float) + block['mrate'] = np.array(mrate_list, dtype=float) + blocks.append(block) + return header, tuple(blocks) - return header, read_blocks(lines) def read_tr(filename): """ read *a.tr file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['MULTIP'], lines = _read_value(lines, int) - block['GAUGE'], lines = _read_value(lines, int) - block['MODE'], lines = _read_value(lines, int) - # read the values - block['lower_index'] = np.zeros(ntrans, dtype=int) - block['lower_2J'] = np.zeros(ntrans, dtype=int) - block['upper_index'] = np.zeros(ntrans, dtype=int) - block['upper_2J'] = np.zeros(ntrans, dtype=int) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['gf'] = np.zeros(ntrans, dtype=float) - block['rate'] = np.zeros(ntrans, dtype=float) - block['multipole'] = np.zeros(ntrans, dtype=float) - block['uta'] = header['uta'] - if len(lines) > 0 and lines[0].strip() != '': - a = lines[0].split() - if len(a) == 10: - block['uta'] = 1 + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['MULTIP'], idx = _read_value(lines, idx, int) + block['GAUGE'], idx = _read_value(lines, idx, int) + block['MODE'], idx = _read_value(lines, idx, int) + uta = header['uta'] + if idx < len(lines) and lines[idx].strip(): + a_check = lines[idx].split() + if len(a_check) == 10: + uta = 1 else: - block['uta'] = 0 - if block['uta'] > 0: - block['sdev'] = np.zeros(ntrans, dtype=float) - block['rci'] = np.zeros(ntrans, dtype=float) - - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks + uta = 0 + block['uta'] = uta + + lower_idx, lower_2j = [], [] + upper_idx, upper_2j = [], [] + delta_e, gf, rate, multipole = [], [], [], [] + sdev, rci = [], [] + + while idx < len(lines): + line = lines[idx] + if not line.strip(): + idx += 1 + break a = line.split() - block['upper_index'][i] = int(a[0]) - block['upper_2J'][i] = int(a[1]) - block['lower_index'][i] = int(a[2]) - block['lower_2J'][i] = int(a[3]) - block['Delta E'][i] = float(a[4]) - if block['uta'] > 0: - block['sdev'][i] = float(a[5]) - block['gf'][i] = float(a[6]) - block['rate'][i] = float(a[7]) - block['multipole'][i] = float(a[8]) - block['rci'][i] = float(a[9]) + upper_idx.append(int(a[0])) + upper_2j.append(int(a[1])) + lower_idx.append(int(a[2])) + lower_2j.append(int(a[3])) + delta_e.append(float(a[4])) + if uta > 0: + sdev.append(float(a[5])) + gf.append(float(a[6])) + rate.append(float(a[7])) + multipole.append(float(a[8])) + rci.append(float(a[9])) else: - block['gf'][i] = float(a[5]) - block['rate'][i] = float(a[6]) - block['multipole'][i] = float(a[7]) - - return (block, ) + gf.append(float(a[5])) + rate.append(float(a[6])) + multipole.append(float(a[7])) + idx += 1 + + block['lower_index'] = np.array(lower_idx, dtype=int) + block['lower_2J'] = np.array(lower_2j, dtype=int) + block['upper_index'] = np.array(upper_idx, dtype=int) + block['upper_2J'] = np.array(upper_2j, dtype=int) + block['Delta E'] = np.array(delta_e, dtype=float) + block['gf'] = np.array(gf, dtype=float) + block['rate'] = np.array(rate, dtype=float) + block['multipole'] = np.array(multipole, dtype=float) + if uta > 0: + block['sdev'] = np.array(sdev, dtype=float) + block['rci'] = np.array(rci, dtype=float) + blocks.append(block) + return header, tuple(blocks) - return header, read_blocks(lines) def read_ai(filename): - """ read *a.ai file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['EMIN'], lines = _read_value(lines, float) - negrid, lines = _read_value(lines, int) - block['EGRID'] = np.zeros(negrid, dtype=float) - for i in range(negrid): - block['EGRID'][i] = float(lines.pop(0)) - # read the values - block['bound_index'] = np.zeros(ntrans, dtype=int) - block['bound_2J'] = np.zeros(ntrans, dtype=int) - block['free_index'] = np.zeros(ntrans, dtype=int) - block['free_2J'] = np.zeros(ntrans, dtype=int) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['AI rate'] = np.zeros(ntrans, dtype=float) - block['DC strength'] = np.zeros(ntrans, dtype=float) - - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['EMIN'], idx = _read_value(lines, idx, float) + negrid, idx = _read_value(lines, idx, int) + block['EGRID'] = np.array([float(lines[i]) for i in range(idx, idx + negrid)], dtype=float) + idx += negrid + + bound_idx, bound_2j, free_idx, free_2j = [], [], [], [] + delta_e, ai_rate, dc_strength = [], [], [] + + while idx < len(lines): + line = lines[idx] + if not line.strip(): + idx += 1 + break a = line.split() - block['bound_index'][i] = int(a[0]) - block['bound_2J'][i] = int(a[1]) - block['free_index'][i] = int(a[2]) - block['free_2J'][i] = int(a[3]) - block['Delta E'][i] = float(a[4]) - block['AI rate'][i] = float(a[5]) - block['DC strength'][i] = float(a[6]) - - return (block, ) - - return header, read_blocks(lines) + bound_idx.append(int(a[0])) + bound_2j.append(int(a[1])) + free_idx.append(int(a[2])) + free_2j.append(int(a[3])) + delta_e.append(float(a[4])) + ai_rate.append(float(a[5])) + dc_strength.append(float(a[6])) + idx += 1 + + block['bound_index'] = np.array(bound_idx, dtype=int) + block['bound_2J'] = np.array(bound_2j, dtype=int) + block['free_index'] = np.array(free_idx, dtype=int) + block['free_2J'] = np.array(free_2j, dtype=int) + block['Delta E'] = np.array(delta_e, dtype=float) + block['AI rate'] = np.array(ai_rate, dtype=float) + block['DC strength'] = np.array(dc_strength, dtype=float) + blocks.append(block) + return header, tuple(blocks) def read_ce(filename): - """ read *a.ce file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['QKMODE'], lines = _read_value(lines, int) - nparams, lines = _read_value(lines, int) - block['MSUB'], lines = _read_value(lines, int) - block['PWTYPE'], lines = _read_value(lines, int) - ntegrid, lines = _read_value(lines, int) - block['TEGRID'] = np.zeros(ntegrid, dtype=float) - for i in range(ntegrid): - block['TEGRID'][i] = float(lines.pop(0)) - - block['TE0'], lines = _read_value(lines, float) - block['ETYPE'], lines = _read_value(lines, int) - negrid, lines = _read_value(lines, int) - block['EGRID'] = np.zeros(negrid, dtype=float) - for i in range(negrid): - block['EGRID'][i] = float(lines.pop(0)) - block['UTYPE'], lines = _read_value(lines, int) - nusr, lines = _read_value(lines, int) - block['USR'] = np.zeros(nusr, dtype=float) - for i in range(nusr): - block['USR'][i] = float(lines.pop(0)) - - block['lower_index'] = np.zeros(ntrans, dtype=int) - block['lower_2J'] = np.zeros(ntrans, dtype=int) - block['upper_index'] = np.zeros(ntrans, dtype=int) - block['upper_2J'] = np.zeros(ntrans, dtype=int) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['bethe'] = np.zeros(ntrans, dtype=float) - block['born'] = np.zeros((ntrans, 2), dtype=float) - if block['MSUB']: - block['collision strength'] = [None] * ntrans - block['crosssection'] = [None] * ntrans - else: - block['collision strength'] = np.zeros((ntrans, nusr), dtype=float) - block['crosssection'] = np.zeros((ntrans, nusr), dtype=float) + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['QKMODE'], idx = _read_value(lines, idx, int) + nparams, idx = _read_value(lines, idx, int) + block['MSUB'], idx = _read_value(lines, idx, int) + block['PWTYPE'], idx = _read_value(lines, idx, int) + ntegrid, idx = _read_value(lines, idx, int) + block['TEGRID'] = np.array([float(lines[i]) for i in range(idx, idx + ntegrid)], dtype=float) + idx += ntegrid + block['TE0'], idx = _read_value(lines, idx, float) + block['ETYPE'], idx = _read_value(lines, idx, int) + negrid, idx = _read_value(lines, idx, int) + block['EGRID'] = np.array([float(lines[i]) for i in range(idx, idx + negrid)], dtype=float) + idx += negrid + block['UTYPE'], idx = _read_value(lines, idx, int) + nusr, idx = _read_value(lines, idx, int) + block['USR'] = np.array([float(lines[i]) for i in range(idx, idx + nusr)], dtype=float) + idx += nusr + lower_idx, lower_2j, upper_idx, upper_2j = [], [], [], [] + delta_e, bethe = [], [] + born_0, born_1 = [], [] if block['MSUB']: - block['ratio collision strength'] = np.zeros(ntrans, dtype=float) + collision_str_all = [] + crosssection_all = [] + ratio_cs = [] + else: + collision_str_all = [] + crosssection_all = [] - nsub = np.zeros(ntrans, dtype=int) - if block['QKMODE'] == 2: - block['params'] == np.zeros((ntrans, 4), dtype=float) + params_all = [] if block['QKMODE'] == 2 else None for tr in range(ntrans): - line = lines[0] - lines = lines[1:] - a = line.split() - block['lower_index'][tr] = int(a[0]) - block['lower_2J'][tr] = int(a[1]) - block['upper_index'][tr] = int(a[2]) - block['upper_2J'][tr] = int(a[3]) - block['Delta E'][tr] = float(a[4]) + a = lines[idx].split() + idx += 1 + lower_idx.append(int(a[0])) + lower_2j.append(int(a[1])) + upper_idx.append(int(a[2])) + upper_2j.append(int(a[3])) + delta_e.append(float(a[4])) nsub = int(a[5]) - if block['MSUB']: - block['collision strength'][tr] = np.zeros( - (nusr, nsub), dtype=float) - block['crosssection'][tr] = np.zeros( - (nusr, nsub), dtype=float) - line = lines[0] - lines = lines[1:] - a = line.split() - block['bethe'][tr] = float(a[0]) - block['born'][tr, 0] = float(a[1]) - block['born'][tr, 1] = float(a[2]) - - for sub in range(nsub): - if block['MSUB']: - line = lines[0] - lines = lines[1:] - block['ratio collision strength'][tr] = float(line) - for i in range(nusr): - line = lines[0] - lines = lines[1:] - a = line.split() - block['collision strength'][tr][i, sub] = float(a[1]) - block['crosssection'][tr][i, sub] = float(a[2]) - if sub < nsub - 1: - line = lines[0] - lines = lines[1:] # skip separator ----- + a_bethe = lines[idx].split() + idx += 1 + bethe.append(float(a_bethe[0])) + born_0.append(float(a_bethe[1])) + born_1.append(float(a_bethe[2])) - else: + if block['QKMODE'] == 2: + params_all.append([float(l) for l in lines[idx].split()]) + idx += 1 + + if block['MSUB']: + tr_ratio_cs = float(lines[idx]) + idx += 1 + ratio_cs.append(tr_ratio_cs) + + tr_cs = [] + tr_xs = [] + for sub in range(nsub): + sub_cs = [] + sub_xs = [] for i in range(nusr): - line = lines[0] - lines = lines[1:] - a = line.split() - block['collision strength'][tr, i] = float(a[1]) - block['crosssection'][tr, i] = float(a[2]) + a_sub = lines[idx].split() + idx += 1 + sub_cs.append(float(a_sub[1])) + sub_xs.append(float(a_sub[2])) + tr_cs.append(sub_cs) + tr_xs.append(sub_xs) + if sub < nsub - 1: + idx += 1 + collision_str_all.append(np.array(tr_cs, dtype=float)) + crosssection_all.append(np.array(tr_xs, dtype=float)) + else: + tr_cs = [] + tr_xs = [] + for i in range(nusr): + a_sub = lines[idx].split() + idx += 1 + tr_cs.append(float(a_sub[1])) + tr_xs.append(float(a_sub[2])) + collision_str_all.append(tr_cs) + crosssection_all.append(tr_xs) + + block['lower_index'] = np.array(lower_idx, dtype=int) + block['lower_2J'] = np.array(lower_2j, dtype=int) + block['upper_index'] = np.array(upper_idx, dtype=int) + block['upper_2J'] = np.array(upper_2j, dtype=int) + block['Delta E'] = np.array(delta_e, dtype=float) + block['bethe'] = np.array(bethe, dtype=float) + block['born'] = np.column_stack((born_0, born_1)) - if len(lines) < 3: - return (block, ) + if block['MSUB']: + block['collision strength'] = collision_str_all + block['crosssection'] = crosssection_all + block['ratio collision strength'] = np.array(ratio_cs, dtype=float) + else: + block['collision strength'] = np.array(collision_str_all, dtype=float) + block['crosssection'] = np.array(crosssection_all, dtype=float) - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks + if block['QKMODE'] == 2: + block['params'] = np.array(params_all, dtype=float) - raise ValueError('Bad file format.') + blocks.append(block) + return header, tuple(blocks) - return header, read_blocks(lines) def read_ci(filename): @@ -817,203 +864,203 @@ def read_ci(filename): with open(filename, 'r') as f: lines = f.readlines() - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] + return header, () + idx += 1 - def read_blocks(lines): + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['QKMODE'], lines = _read_value(lines, int) - nparams, lines = _read_value(lines, int) - block['PWTYPE'], lines = _read_value(lines, int) - ntegrid, lines = _read_value(lines, int) - block['TEGRID'] = np.zeros(ntegrid, dtype=float) - for i in range(ntegrid): - block['TEGRID'][i] = float(lines.pop(0)) - block['ETYPE'], lines = _read_value(lines, int) - negrid, lines = _read_value(lines, int) - block['EGRID'] = np.zeros(negrid, dtype=float) - for i in range(negrid): - block['EGRID'][i] = float(lines.pop(0)) - block['UTYPE'], lines = _read_value(lines, int) - nusr, lines = _read_value(lines, int) - block['USR'] = np.zeros(nusr, dtype=float) - for i in range(nusr): - block['USR'][i] = float(lines.pop(0)) - - # read the values - block['bound_index'] = np.zeros(ntrans, dtype=int) - block['bound_2J'] = np.zeros(ntrans, dtype=int) - block['free_index'] = np.zeros(ntrans, dtype=int) - block['free_2J'] = np.zeros(ntrans, dtype=int) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['Delta L'] = np.zeros(ntrans, dtype=int) - block['parameters'] = np.zeros((ntrans, nparams), dtype=float) - block['collision strength'] = np.zeros((ntrans, nusr), dtype=float) - block['crosssection'] = np.zeros((ntrans, nusr), dtype=float) + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['QKMODE'], idx = _read_value(lines, idx, int) + nparams, idx = _read_value(lines, idx, int) + block['PWTYPE'], idx = _read_value(lines, idx, int) + ntegrid, idx = _read_value(lines, idx, int) + block['TEGRID'] = np.array([float(lines[i]) for i in range(idx, idx + ntegrid)], dtype=float) + idx += ntegrid + block['ETYPE'], idx = _read_value(lines, idx, int) + negrid, idx = _read_value(lines, idx, int) + block['EGRID'] = np.array([float(lines[i]) for i in range(idx, idx + negrid)], dtype=float) + idx += negrid + block['UTYPE'], idx = _read_value(lines, idx, int) + nusr, idx = _read_value(lines, idx, int) + block['USR'] = np.array([float(lines[i]) for i in range(idx, idx + nusr)], dtype=float) + idx += nusr + + bound_idx, bound_2j, free_idx, free_2j = [], [], [], [] + delta_e, delta_l = [], [] + parameters_list = [] + collision_str_list = [] + crosssection_list = [] for tr in range(ntrans): - line = lines[0] - lines = lines[1:] - a = line.split() - block['bound_index'][tr] = int(a[0]) - block['bound_2J'][tr] = int(a[1]) - block['free_index'][tr] = int(a[2]) - block['free_2J'][tr] = int(a[3]) - block['Delta E'][tr] = float(a[4]) - block['Delta L'][tr] = int(a[5]) - block['parameters'][tr] = [float(l) for l in lines[0].split()] - lines = lines[1:] + a = lines[idx].split() + idx += 1 + bound_idx.append(int(a[0])) + bound_2j.append(int(a[1])) + free_idx.append(int(a[2])) + free_2j.append(int(a[3])) + delta_e.append(float(a[4])) + delta_l.append(int(a[5])) + param_line = lines[idx].split() + idx += 1 + parameters_list.append([float(l) for l in param_line[:nparams]]) + cs_row = [] + xs_row = [] for i in range(nusr): - line = lines[0] - lines = lines[1:] - a = line.split() - block['collision strength'][tr, i] = float(a[1]) - block['crosssection'][tr, i] = float(a[2]) + a_sub = lines[idx].split() + idx += 1 + # a_sub[0] is the collision energy (egrid) + cs_row.append(float(a_sub[1])) + xs_row.append(float(a_sub[2])) + collision_str_list.append(cs_row) + crosssection_list.append(xs_row) - if len(lines) < 3: - return (block, ) + block['bound_index'] = np.array(bound_idx, dtype=int) + block['bound_2J'] = np.array(bound_2j, dtype=int) + block['free_index'] = np.array(free_idx, dtype=int) + block['free_2J'] = np.array(free_2j, dtype=int) + block['Delta E'] = np.array(delta_e, dtype=float) + block['Delta L'] = np.array(delta_l, dtype=int) + block['parameters'] = np.array(parameters_list, dtype=float) + block['collision strength'] = np.array(collision_str_list, dtype=float) + block['crosssection'] = np.array(crosssection_list, dtype=float) - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks + blocks.append(block) - raise ValueError('Bad file format.') - - return header, read_blocks(lines) + return header, tuple(blocks) def read_rr(filename): - """ read *a.rr file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['QKMODE'], lines = _read_value(lines, int) - block['MULTIP'], lines = _read_value(lines, int) - nparams, lines = _read_value(lines, int) - ntegrid, lines = _read_value(lines, int) - block['TEGRID'] = np.zeros(ntegrid, dtype=float) - for i in range(ntegrid): - block['TEGRID'][i] = float(lines.pop(0)) - block['ETYPE'], lines = _read_value(lines, int) - negrid, lines = _read_value(lines, int) - block['EGRID'] = np.zeros(negrid, dtype=float) - for i in range(negrid): - block['EGRID'][i] = float(lines.pop(0)) - block['UTYPE'], lines = _read_value(lines, int) - nusr, lines = _read_value(lines, int) - block['USR'] = np.zeros(nusr, dtype=float) - for i in range(nusr): - block['USR'][i] = float(lines.pop(0)) - - # read the values - block['bound_index'] = np.zeros(ntrans, dtype=int) - block['bound_2J'] = np.zeros(ntrans, dtype=int) - block['free_index'] = np.zeros(ntrans, dtype=int) - block['free_2J'] = np.zeros(ntrans, dtype=int) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['Delta L'] = np.zeros(ntrans, dtype=int) - block['parameters'] = np.zeros((ntrans, nparams), dtype=float) - block['RR crosssection'] = np.zeros((ntrans, nusr), dtype=float) - block['PI crosssection'] = np.zeros((ntrans, nusr), dtype=float) - block['gf'] = np.zeros((ntrans, nusr), dtype=float) + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['QKMODE'], idx = _read_value(lines, idx, int) + block['MULTIP'], idx = _read_value(lines, idx, int) + nparams, idx = _read_value(lines, idx, int) + + ntegrid, idx = _read_value(lines, idx, int) + block['TEGRID'] = np.array([float(lines[i]) for i in range(idx, idx + ntegrid)], dtype=float) + idx += ntegrid + + block['ETYPE'], idx = _read_value(lines, idx, int) + negrid, idx = _read_value(lines, idx, int) + block['EGRID'] = np.array([float(lines[i]) for i in range(idx, idx + negrid)], dtype=float) + idx += negrid + + block['UTYPE'], idx = _read_value(lines, idx, int) + nusr, idx = _read_value(lines, idx, int) + block['USR'] = np.array([float(lines[i]) for i in range(idx, idx + nusr)], dtype=float) + idx += nusr + + bound_idx, bound_2j, free_idx, free_2j = [], [], [], [] + delta_e, delta_l = [], [] + parameters_list = [] + rr_cs_list = [] + pi_cs_list = [] + gf_list = [] for tr in range(ntrans): - line = lines[0] - lines = lines[1:] - a = line.split() - block['bound_index'][tr] = int(a[0]) - block['bound_2J'][tr] = int(a[1]) - block['free_index'][tr] = int(a[2]) - block['free_2J'][tr] = int(a[3]) - block['Delta E'][tr] = float(a[4]) - block['Delta L'][tr] = int(a[5]) - block['parameters'][tr] = [float(l) for l in lines[0].split()] - lines = lines[1:] + a = lines[idx].split() + idx += 1 + bound_idx.append(int(a[0])) + bound_2j.append(int(a[1])) + free_idx.append(int(a[2])) + free_2j.append(int(a[3])) + delta_e.append(float(a[4])) + delta_l.append(int(a[5])) + + parameters_list.append([float(l) for l in lines[idx].split()]) + idx += 1 + + tr_rr = [] + tr_pi = [] + tr_gf = [] for i in range(nusr): - line = lines[0] - lines = lines[1:] - a = line.split() - block['RR crosssection'][tr, i] = float(a[1]) - block['PI crosssection'][tr, i] = float(a[2]) - block['gf'][tr, i] = float(a[3]) - - if len(lines) < 3: - return (block, ) - - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks - - raise ValueError('Bad file format.') - - return header, read_blocks(lines) + a_sub = lines[idx].split() + idx += 1 + tr_rr.append(float(a_sub[1])) + tr_pi.append(float(a_sub[2])) + tr_gf.append(float(a_sub[3])) + + rr_cs_list.append(tr_rr) + pi_cs_list.append(tr_pi) + gf_list.append(tr_gf) + + block['bound_index'] = np.array(bound_idx, dtype=int) + block['bound_2J'] = np.array(bound_2j, dtype=int) + block['free_index'] = np.array(free_idx, dtype=int) + block['free_2J'] = np.array(free_2j, dtype=int) + block['Delta E'] = np.array(delta_e, dtype=float) + block['Delta L'] = np.array(delta_l, dtype=int) + block['parameters'] = np.array(parameters_list, dtype=float) + block['RR crosssection'] = np.array(rr_cs_list, dtype=float) + block['PI crosssection'] = np.array(pi_cs_list, dtype=float) + block['gf'] = np.array(gf_list, dtype=float) + blocks.append(block) + return header, tuple(blocks) def read_sp(filename): - """ read *a.sp file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - block['NELE'], lines = _read_value(lines, int) - ntrans, lines = _read_value(lines, int) - block['TYPE'], lines = _read_value(lines, str) - block['IBLK'], lines = _read_value(lines, int) - block['ICOMP'], lines = _read_value(lines, str) - block['FBLK'], lines = _read_value(lines, int) - block['FCOMP'], lines = _read_value(lines, str) - - # read the values - block['block'] = np.zeros(ntrans, dtype=int) - block['level'] = np.zeros(ntrans, dtype=int) - block['abs. energy'] = np.zeros(ntrans, dtype=float) - block['population'] = np.zeros(ntrans, dtype=float) - block['Delta E'] = np.zeros(ntrans, dtype=float) - block['emissivity'] = np.zeros(ntrans, dtype=float) + block['NELE'], idx = _read_value(lines, idx, int) + ntrans, idx = _read_value(lines, idx, int) + block['TYPE'], idx = _read_value(lines, idx, str) + block['IBLK'], idx = _read_value(lines, idx, int) + block['ICOMP'], idx = _read_value(lines, idx, str) + block['FBLK'], idx = _read_value(lines, idx, int) + block['FCOMP'], idx = _read_value(lines, idx, str) + + blk, level = [], [] + abs_energy, population = [], [] + delta_e, emissivity = [], [] for tr in range(ntrans): - line = lines[0] - lines = lines[1:] - a = line.split() - block['block'][tr] = int(a[0]) - block['level'][tr] = int(a[1]) - block['abs. energy'][tr] = float(a[2]) - block['population'][tr] = float(a[3]) - block['Delta E'][tr] = float(a[4]) - block['emissivity'][tr] = float(a[5]) - - for i, line in enumerate(lines): - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks - - return (block, ) + a = lines[idx].split() + idx += 1 + blk.append(int(a[0])) + level.append(int(a[1])) + abs_energy.append(float(a[2])) + population.append(float(a[3])) + delta_e.append(float(a[4])) + emissivity.append(float(a[5])) + + block['block'] = np.array(blk, dtype=int) + block['level'] = np.array(level, dtype=int) + block['abs. energy'] = np.array(abs_energy, dtype=float) + block['population'] = np.array(population, dtype=float) + block['Delta E'] = np.array(delta_e, dtype=float) + block['emissivity'] = np.array(emissivity, dtype=float) + blocks.append(block) + return header, tuple(blocks) - return header, read_blocks(lines) def read_wfun(fn, npi=0, rmax=None): r = np.loadtxt(fn, unpack=1) @@ -1116,7 +1163,7 @@ def read_den(fn, cfg=1, rnd=4): else: nq = 0.0 return cfgnr(nlq) - + def read_pot(fn, cfg=None, header=None, rnd=4): if cfg is None and header is None: return np.loadtxt(fn, unpack=1) @@ -1341,72 +1388,66 @@ def interp_rra(d, ea, aa=None): r[3,i] = b1 r[4,i] = -b1/b0 return r - + + def read_rt(filename): - """ read *a.rt file. """ with open(filename, 'r') as f: lines = f.readlines() - - # header - header, lines = _get_header(lines) + header, idx = _get_header(lines, idx=0) if header['NBlocks'] == 0: - return header,() - lines = lines[1:] - - def read_blocks(lines): + return header, () + idx += 1 + blocks = [] + while idx < len(lines): + if not lines[idx].strip(): + idx += 1 + continue block = {} - ntrans, lines = _read_value(lines, int) - block['EDEN'], lines = _read_value(lines, float) - block['EDIST'], lines = _read_value(lines, int) - npedis, lines = _read_value(lines, int) - block['EDIS'] = np.zeros(npedis, float) - for i in range(npedis): - line = lines[0] - lines = lines[1:] - block['EDIS'][i] = float(line) - block['PDEN'], lines = _read_value(lines, float) - block['PDIST'], lines = _read_value(lines, int) - nppdis, lines = _read_value(lines, int) - block['PPDIS'] = np.zeros(nppdis, float) - for i in range(npedis): - line = lines[0] - lines = lines[1:] - block['PPDIS'][i] = float(line) - lines = lines[1:] # skip header - - # read the values - block['block'] = np.zeros(ntrans, dtype=int) - block['level'] = np.zeros(ntrans, dtype=int) - block['NB'] = np.zeros(ntrans, dtype=int) - block['TR'] = np.zeros(ntrans, dtype=int) - block['CE'] = np.zeros(ntrans, dtype=int) - block['RR'] = np.zeros(ntrans, dtype=int) - block['AI'] = np.zeros(ntrans, dtype=int) - block['CI'] = np.zeros(ntrans, dtype=int) - block['ncomplex'] = np.chararray(ntrans, itemsize=20) - - for tr in range(ntrans): - line = lines[0] - lines = lines[1:] - - if line.strip() == '': # if empty - blocks = read_blocks(lines[i+1:]) - return (block, ) + blocks - - block['block'][tr] = int(line[:6]) - block['level'][tr] = int(line[7:11]) - block['NB'][tr] = float(line[12:24]) - block['TR'][tr] = float(line[25:36]) - block['CE'][tr] = float(line[37:48]) - block['RR'][tr] = float(line[49:60]) - block['AI'][tr] = float(line[61:72]) - block['CI'][tr] = float(line[73:84]) - block['ncomplex'][tr] = line[85:].strip() - - return (block, ) - - return header, read_blocks(lines) - + ntrans, idx = _read_value(lines, idx, int) + block['EDEN'], idx = _read_value(lines, idx, float) + block['EDIST'], idx = _read_value(lines, idx, int) + npedis, idx = _read_value(lines, idx, int) + block['EDIS'] = np.array([float(lines[i]) for i in range(idx, idx + npedis)], dtype=float) + idx += npedis + + block['PDEN'], idx = _read_value(lines, idx, float) + block['PDIST'], idx = _read_value(lines, idx, int) + nppdis, idx = _read_value(lines, idx, int) + block['PPDIS'] = np.array([float(lines[i]) for i in range(idx, idx + nppdis)], dtype=float) + idx += nppdis + + idx += 1 + + blk, level = [], [] + nb, tr, ce, rr, ai, ci = [], [], [], [], [], [] + ncomplex = [] + + for tr_idx in range(ntrans): + line = lines[idx] + idx += 1 + if not line.strip(): + break + blk.append(int(line[:6])) + level.append(int(line[7:11])) + nb.append(float(line[12:24])) + tr.append(float(line[25:36])) + ce.append(float(line[37:48])) + rr.append(float(line[49:60])) + ai.append(float(line[61:72])) + ci.append(float(line[73:84])) + ncomplex.append(line[85:].strip()) + + block['block'] = np.array(blk, dtype=int) + block['level'] = np.array(level, dtype=int) + block['NB'] = np.array(nb, dtype=float) + block['TR'] = np.array(tr, dtype=float) + block['CE'] = np.array(ce, dtype=float) + block['RR'] = np.array(rr, dtype=float) + block['AI'] = np.array(ai, dtype=float) + block['CI'] = np.array(ci, dtype=float) + block['ncomplex'] = np.array(ncomplex, dtype='U20') + blocks.append(block) + return header, tuple(blocks) MAX_SYMMETRIES = 256 From e72c5c04807ab923880f338cabe4655f05993f70 Mon Sep 17 00:00:00 2001 From: Mingyu Zhu Date: Tue, 2 Jun 2026 20:03:22 +0800 Subject: [PATCH 2/4] fixed `rfac.read_ce` function when handling `MSUB=1`. --- python/rfac.py | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/python/rfac.py b/python/rfac.py index e7d665a7..67a23b70 100644 --- a/python/rfac.py +++ b/python/rfac.py @@ -769,18 +769,14 @@ def read_ce(filename): nusr, idx = _read_value(lines, idx, int) block['USR'] = np.array([float(lines[i]) for i in range(idx, idx + nusr)], dtype=float) idx += nusr + lower_idx, lower_2j, upper_idx, upper_2j = [], [], [], [] delta_e, bethe = [], [] born_0, born_1 = [], [] - if block['MSUB']: - collision_str_all = [] - crosssection_all = [] - ratio_cs = [] - else: - collision_str_all = [] - crosssection_all = [] - + collision_str_all = [] + crosssection_all = [] + ratio_cs = [] params_all = [] if block['QKMODE'] == 2 else None for tr in range(ntrans): @@ -804,15 +800,17 @@ def read_ce(filename): idx += 1 if block['MSUB']: - tr_ratio_cs = float(lines[idx]) - idx += 1 - ratio_cs.append(tr_ratio_cs) - + tr_ratio_cs = [] tr_cs = [] tr_xs = [] for sub in range(nsub): + + tr_ratio_cs.append(float(lines[idx].strip())) + idx += 1 + sub_cs = [] sub_xs = [] + for i in range(nusr): a_sub = lines[idx].split() idx += 1 @@ -820,10 +818,13 @@ def read_ce(filename): sub_xs.append(float(a_sub[2])) tr_cs.append(sub_cs) tr_xs.append(sub_xs) + if sub < nsub - 1: idx += 1 - collision_str_all.append(np.array(tr_cs, dtype=float)) - crosssection_all.append(np.array(tr_xs, dtype=float)) + + collision_str_all.append(np.array(tr_cs, dtype=float).T) + crosssection_all.append(np.array(tr_xs, dtype=float).T) + ratio_cs.append(tr_ratio_cs) else: tr_cs = [] tr_xs = [] @@ -846,6 +847,7 @@ def read_ce(filename): if block['MSUB']: block['collision strength'] = collision_str_all block['crosssection'] = crosssection_all + # check this because in the old code only the first value of each group is picked block['ratio collision strength'] = np.array(ratio_cs, dtype=float) else: block['collision strength'] = np.array(collision_str_all, dtype=float) From 11f2e6d21d0a263db4d5ffe63d7b76dccc4b6e20 Mon Sep 17 00:00:00 2001 From: Mingyu Zhu Date: Wed, 2 Sep 2026 14:04:36 +0800 Subject: [PATCH 3/4] fix faclib/makefile in f2c compilation flag. --- faclib/Makefile.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/faclib/Makefile.in b/faclib/Makefile.in index 3c1c7870..ec8488af 100644 --- a/faclib/Makefile.in +++ b/faclib/Makefile.in @@ -26,7 +26,7 @@ endif faclib: ${OBJS} ifdef F2C - $(FC) -c -o f2c.o f2c.f90 + $(FC) $(FFLAGS) -c -o f2c.o f2c.f90 endif ar r ${TOPDIR}/libfac.a ${OBJS} ${FOBJ} From 46b622624441e7371b750c93a5e5672ab7e35636 Mon Sep 17 00:00:00 2001 From: Mingyu Zhu Date: Fri, 4 Sep 2026 21:30:05 +0800 Subject: [PATCH 4/4] fix int overflow in hamiltonian matrix size calculation in structure.c --- faclib/structure.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/faclib/structure.c b/faclib/structure.c index 5a428b82..c87e7255 100644 --- a/faclib/structure.c +++ b/faclib/structure.c @@ -460,7 +460,7 @@ int ConstructHamiltonDiagonal(int isym, int k, int *kg, int m) { } if (!(h->basis)) goto ERROR; - h->msize = h->dim*h->n_basis + h->dim; + h->msize = (size_t)h->dim * (size_t)h->n_basis + h->dim; if (h->hamilton == NULL) { h->hsize0 = h->hsize; h->hamilton = (double *) malloc(sizeof(double)*h->hsize); @@ -7090,7 +7090,7 @@ int AllocHamMem(HAMILTON *h, int hdim, int nbasis) { h->iham = nhams; h->n_basis = nbasis; - h->msize = h->dim * h->n_basis + h->dim; + h->msize = (size_t)h->dim * (size_t)h->n_basis + (size_t)h->dim; if (h->mixing == NULL) { h->msize0 = h->msize; h->mixing = (double *) malloc(sizeof(double)*(size_t)h->msize); @@ -7127,10 +7127,10 @@ int AllocHamMem(HAMILTON *h, int hdim, int nbasis) { return 0; } - t = hdim*(hdim+1)/2; + t = (size_t)hdim * (size_t)(hdim+1)/2; h->dsize = t; h->dsize2 = t*2; - h->hsize = t + hdim*jp + jp; + h->hsize = t + (size_t)hdim*(size_t)jp + jp; if (h->hamilton == NULL) { h->hsize0 = h->hsize; h->hamilton = (double *) malloc(sizeof(double)*(size_t)h->hsize); @@ -7147,7 +7147,7 @@ int AllocHamMem(HAMILTON *h, int hdim, int nbasis) { size_t wl = 3*t; size_t wi = 0; if (nbasis > hdim) { - wl += t + hdim*(nbasis-hdim); + wl += t + (size_t)hdim*(size_t)(nbasis-hdim); wi += nbasis-hdim; } if (h->work == NULL) {