scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-10-05 07:45:59 +02:00
Code Issues Releases Wiki Activity

Allow no basis set

Browse Source
This commit is contained in:
Anthony Scemama 2023-05-10 14:44:30 +02:00
parent 192f4c6913
commit 6f04cc0eb0
2 changed files with 99 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
devel/trexio/import_trexio_integrals.irp.f
View File

@ -1,8 +1,22 @@
program import_integrals_ao 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 end
subroutine run subroutine run(f)
use trexio use trexio
use map_module use map_module
implicit none implicit none
@ -10,7 +24,7 @@ subroutine run
! Program to import integrals from TREXIO ! Program to import integrals from TREXIO
END_DOC END_DOC
integer(trexio_t) :: f ! TREXIO file handle integer(trexio_t), intent(in) :: f ! TREXIO file handle
integer(trexio_exit_code) :: rc integer(trexio_exit_code) :: rc
integer ::i,j,k,l integer ::i,j,k,l
@ -25,16 +39,6 @@ subroutine run
double precision, allocatable :: V(:) double precision, allocatable :: V(:)
integer , allocatable :: Vi(:,:) integer , allocatable :: Vi(:,:)
double precision :: s 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 if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
rc = trexio_read_nucleus_repulsion(f, s) rc = trexio_read_nucleus_repulsion(f, s)

147
devel/trexio/qp_import_trexio.py
View File

@ -72,22 +72,30 @@ def write_ezfio(trexio_filename, filename):
print("Nuclei\t\t...\t", end=' ') print("Nuclei\t\t...\t", end=' ')
charge = trexio.read_nucleus_charge(trexio_file) charge = [0.]
ezfio.set_nuclei_nucl_num(len(charge)) if trexio.has_nucleus(trexio_file):
ezfio.set_nuclei_nucl_charge(charge) 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 = trexio.read_nucleus_coord(trexio_file)
coord = np.transpose(coord) coord = np.transpose(coord)
ezfio.set_nuclei_nucl_coord(coord) ezfio.set_nuclei_nucl_coord(coord)
label = trexio.read_nucleus_label(trexio_file) label = trexio.read_nucleus_label(trexio_file)
nucl_num = trexio.read_nucleus_num(trexio_file) nucl_num = trexio.read_nucleus_num(trexio_file)
# Transformt H1 into H # Transformt H1 into H
import re import re
p = re.compile(r'(\d*)$') p = re.compile(r'(\d*)$')
label = [p.sub("", x).capitalize() for x in label] label = [p.sub("", x).capitalize() for x in label]
ezfio.set_nuclei_nucl_label(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") print("OK")
@ -104,6 +112,9 @@ def write_ezfio(trexio_filename, filename):
except: except:
num_alpha = sum(charge) - num_beta 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_alpha_num(num_alpha)
ezfio.set_electrons_elec_beta_num(num_beta) ezfio.set_electrons_elec_beta_num(num_beta)
@ -111,10 +122,11 @@ def write_ezfio(trexio_filename, filename):
print("Basis\t\t...\t", end=' ') print("Basis\t\t...\t", end=' ')
shell_num = 0
try: try:
basis_type = trexio.read_basis_type(trexio_file) basis_type = trexio.read_basis_type(trexio_file)
if basis_type.lower() != "gaussian": if basis_type.lower() not in ["gaussian", "slater"]:
raise TypeError raise TypeError
shell_num = trexio.read_basis_shell_num(trexio_file) 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=' ') print("AOS\t\t...\t", end=' ')
ao_shell = trexio.read_ao_shell(trexio_file) try:
cartesian = trexio.read_ao_cartesian(trexio_file) cartesian = trexio.read_ao_cartesian(trexio_file)
except:
cartesian = True
if not cartesian: if not cartesian:
raise TypeError('Only cartesian TREXIO files can be converted') raise TypeError('Only cartesian TREXIO files can be converted')
ao_num = trexio.read_ao_num(trexio_file) ao_num = trexio.read_ao_num(trexio_file)
ezfio.set_ao_basis_ao_num(ao_num) ezfio.set_ao_basis_ao_num(ao_num)
at = [ nucl_index[i]+1 for i in ao_shell ] if shell_num > 0:
ezfio.set_ao_basis_ao_nucl(at) 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) ] num_prim0 = [ 0 for i in range(shell_num) ]
for i in shell_index: for i in shell_index:
num_prim0[i] += 1 num_prim0[i] += 1
coef = {} coef = {}
expo = {} expo = {}
for i,c in enumerate(coefficient): for i,c in enumerate(coefficient):
idx = shell_index[i] idx = shell_index[i]
if idx in coef: if idx in coef:
coef[idx].append(c) coef[idx].append(c)
expo[idx].append(exponent[i]) expo[idx].append(exponent[i])
else: else:
coef[idx] = [c] coef[idx] = [c]
expo[idx] = [exponent[i]] expo[idx] = [exponent[i]]
coefficient = [] coefficient = []
exponent = [] exponent = []
power_x = [] power_x = []
power_y = [] power_y = []
power_z = [] power_z = []
num_prim = [] num_prim = []
for i in range(shell_num): for i in range(shell_num):
for x,y,z in generate_xyz(ang_mom[i]): for x,y,z in generate_xyz(ang_mom[i]):
power_x.append(x) power_x.append(x)
power_y.append(y) power_y.append(y)
power_z.append(z) power_z.append(z)
coefficient.append(coef[i]) coefficient.append(coef[i])
exponent.append(expo[i]) exponent.append(expo[i])
num_prim.append(num_prim0[i]) num_prim.append(num_prim0[i])
assert (len(coefficient) == ao_num) assert (len(coefficient) == ao_num)
ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
ezfio.set_ao_basis_ao_prim_num(num_prim) 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): for i in range(ao_num):
coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] 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)] exponent [i] += [0. for j in range(len(exponent[i]), prim_num_max)]
coefficient = reduce(lambda x, y: x + y, coefficient, []) coefficient = reduce(lambda x, y: x + y, coefficient, [])
exponent = reduce(lambda x, y: x + y, exponent , []) exponent = reduce(lambda x, y: x + y, exponent , [])
coef = [] coef = []
expo = [] expo = []
for i in range(prim_num_max): for i in range(prim_num_max):
for j in range(i, len(coefficient), prim_num_max): for j in range(i, len(coefficient), prim_num_max):
coef.append(coefficient[j]) coef.append(coefficient[j])
expo.append(exponent[j]) expo.append(exponent[j])
# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max) # ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
ezfio.set_ao_basis_ao_coef(coef) ezfio.set_ao_basis_ao_coef(coef)
ezfio.set_ao_basis_ao_expo(expo) ezfio.set_ao_basis_ao_expo(expo)
ezfio.set_ao_basis_ao_basis("Read from TREXIO") ezfio.set_ao_basis_ao_basis("Read from TREXIO")
print("OK") print("OK")