From 6f04cc0eb0b4e8688bc735eeb93aaaca28b8e4ab Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 10 May 2023 14:44:30 +0200 Subject: [PATCH] Allow no basis set --- devel/trexio/import_trexio_integrals.irp.f | 30 +++-- devel/trexio/qp_import_trexio.py | 147 ++++++++++++--------- 2 files changed, 99 insertions(+), 78 deletions(-) diff --git a/devel/trexio/import_trexio_integrals.irp.f b/devel/trexio/import_trexio_integrals.irp.f index cf610e3..9f9ad9d 100644 --- a/devel/trexio/import_trexio_integrals.irp.f +++ b/devel/trexio/import_trexio_integrals.irp.f @@ -1,8 +1,22 @@ program import_integrals_ao - call run + use trexio + implicit none + integer(trexio_t) :: f ! TREXIO file handle + integer(trexio_exit_code) :: rc + + f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc) + if (f == 0_8) then + print *, 'Unable to open TREXIO file for reading' + print *, 'rc = ', rc + stop -1 + endif + + call run(f) + rc = trexio_close(f) + call trexio_assert(rc, TREXIO_SUCCESS) end -subroutine run +subroutine run(f) use trexio use map_module implicit none @@ -10,7 +24,7 @@ subroutine run ! Program to import integrals from TREXIO END_DOC - integer(trexio_t) :: f ! TREXIO file handle + integer(trexio_t), intent(in) :: f ! TREXIO file handle integer(trexio_exit_code) :: rc integer ::i,j,k,l @@ -25,16 +39,6 @@ subroutine run double precision, allocatable :: V(:) integer , allocatable :: Vi(:,:) double precision :: s - PROVIDE ao_num - - f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc) - if (f == 0_8) then - print *, 'Unable to open TREXIO file for reading' - print *, 'rc = ', rc - stop -1 - endif - - if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then rc = trexio_read_nucleus_repulsion(f, s) diff --git a/devel/trexio/qp_import_trexio.py b/devel/trexio/qp_import_trexio.py index 165fb9d..399ebac 100755 --- a/devel/trexio/qp_import_trexio.py +++ b/devel/trexio/qp_import_trexio.py @@ -72,22 +72,30 @@ def write_ezfio(trexio_filename, filename): print("Nuclei\t\t...\t", end=' ') - charge = trexio.read_nucleus_charge(trexio_file) - ezfio.set_nuclei_nucl_num(len(charge)) - ezfio.set_nuclei_nucl_charge(charge) + charge = [0.] + if trexio.has_nucleus(trexio_file): + charge = trexio.read_nucleus_charge(trexio_file) + ezfio.set_nuclei_nucl_num(len(charge)) + ezfio.set_nuclei_nucl_charge(charge) - coord = trexio.read_nucleus_coord(trexio_file) - coord = np.transpose(coord) - ezfio.set_nuclei_nucl_coord(coord) + coord = trexio.read_nucleus_coord(trexio_file) + coord = np.transpose(coord) + ezfio.set_nuclei_nucl_coord(coord) - label = trexio.read_nucleus_label(trexio_file) - nucl_num = trexio.read_nucleus_num(trexio_file) + label = trexio.read_nucleus_label(trexio_file) + nucl_num = trexio.read_nucleus_num(trexio_file) - # Transformt H1 into H - import re - p = re.compile(r'(\d*)$') - label = [p.sub("", x).capitalize() for x in label] - ezfio.set_nuclei_nucl_label(label) + # Transformt H1 into H + import re + p = re.compile(r'(\d*)$') + label = [p.sub("", x).capitalize() for x in label] + ezfio.set_nuclei_nucl_label(label) + + else: + ezfio.set_nuclei_nucl_num(1) + ezfio.set_nuclei_nucl_charge([0.]) + ezfio.set_nuclei_nucl_coord([0.,0.,0.]) + ezfio.set_nuclei_nucl_label(["X"]) print("OK") @@ -104,6 +112,9 @@ def write_ezfio(trexio_filename, filename): except: num_alpha = sum(charge) - num_beta + if num_alpha == 0: + print("\n\nError: There are zero electrons in the TREXIO file.\n\n") + sys.exit(1) ezfio.set_electrons_elec_alpha_num(num_alpha) ezfio.set_electrons_elec_beta_num(num_beta) @@ -111,10 +122,11 @@ def write_ezfio(trexio_filename, filename): print("Basis\t\t...\t", end=' ') + shell_num = 0 try: basis_type = trexio.read_basis_type(trexio_file) - - if basis_type.lower() != "gaussian": + + if basis_type.lower() not in ["gaussian", "slater"]: raise TypeError shell_num = trexio.read_basis_shell_num(trexio_file) @@ -173,72 +185,77 @@ def write_ezfio(trexio_filename, filename): print("AOS\t\t...\t", end=' ') - ao_shell = trexio.read_ao_shell(trexio_file) - cartesian = trexio.read_ao_cartesian(trexio_file) + try: + cartesian = trexio.read_ao_cartesian(trexio_file) + except: + cartesian = True + if not cartesian: raise TypeError('Only cartesian TREXIO files can be converted') ao_num = trexio.read_ao_num(trexio_file) ezfio.set_ao_basis_ao_num(ao_num) - at = [ nucl_index[i]+1 for i in ao_shell ] - ezfio.set_ao_basis_ao_nucl(at) + if shell_num > 0: + ao_shell = trexio.read_ao_shell(trexio_file) + at = [ nucl_index[i]+1 for i in ao_shell ] + ezfio.set_ao_basis_ao_nucl(at) - num_prim0 = [ 0 for i in range(shell_num) ] - for i in shell_index: - num_prim0[i] += 1 + num_prim0 = [ 0 for i in range(shell_num) ] + for i in shell_index: + num_prim0[i] += 1 - coef = {} - expo = {} - for i,c in enumerate(coefficient): - idx = shell_index[i] - if idx in coef: - coef[idx].append(c) - expo[idx].append(exponent[i]) - else: - coef[idx] = [c] - expo[idx] = [exponent[i]] + coef = {} + expo = {} + for i,c in enumerate(coefficient): + idx = shell_index[i] + if idx in coef: + coef[idx].append(c) + expo[idx].append(exponent[i]) + else: + coef[idx] = [c] + expo[idx] = [exponent[i]] - coefficient = [] - exponent = [] - power_x = [] - power_y = [] - power_z = [] - num_prim = [] + coefficient = [] + exponent = [] + power_x = [] + power_y = [] + power_z = [] + num_prim = [] - for i in range(shell_num): - for x,y,z in generate_xyz(ang_mom[i]): - power_x.append(x) - power_y.append(y) - power_z.append(z) - coefficient.append(coef[i]) - exponent.append(expo[i]) - num_prim.append(num_prim0[i]) + for i in range(shell_num): + for x,y,z in generate_xyz(ang_mom[i]): + power_x.append(x) + power_y.append(y) + power_z.append(z) + coefficient.append(coef[i]) + exponent.append(expo[i]) + num_prim.append(num_prim0[i]) - assert (len(coefficient) == ao_num) - ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) - ezfio.set_ao_basis_ao_prim_num(num_prim) + assert (len(coefficient) == ao_num) + ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) + ezfio.set_ao_basis_ao_prim_num(num_prim) - prim_num_max = max( [ len(x) for x in coefficient ] ) + prim_num_max = max( [ len(x) for x in coefficient ] ) - for i in range(ao_num): - coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] - exponent [i] += [0. for j in range(len(exponent[i]), prim_num_max)] + for i in range(ao_num): + coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] + exponent [i] += [0. for j in range(len(exponent[i]), prim_num_max)] - coefficient = reduce(lambda x, y: x + y, coefficient, []) - exponent = reduce(lambda x, y: x + y, exponent , []) + coefficient = reduce(lambda x, y: x + y, coefficient, []) + exponent = reduce(lambda x, y: x + y, exponent , []) - coef = [] - expo = [] - for i in range(prim_num_max): - for j in range(i, len(coefficient), prim_num_max): - coef.append(coefficient[j]) - expo.append(exponent[j]) + coef = [] + expo = [] + for i in range(prim_num_max): + for j in range(i, len(coefficient), prim_num_max): + coef.append(coefficient[j]) + expo.append(exponent[j]) -# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max) - ezfio.set_ao_basis_ao_coef(coef) - ezfio.set_ao_basis_ao_expo(expo) - ezfio.set_ao_basis_ao_basis("Read from TREXIO") +# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max) + ezfio.set_ao_basis_ao_coef(coef) + ezfio.set_ao_basis_ao_expo(expo) + ezfio.set_ao_basis_ao_basis("Read from TREXIO") print("OK")