Added TREXIO module

This commit is contained in:
Anthony Scemama 2023-05-11 22:48:48 +02:00
parent 01b70ffb17
commit 4959882293
9 changed files with 1332 additions and 0 deletions

54
src/trexio/EZFIO.cfg Normal file
View File

@ -0,0 +1,54 @@
[backend]
type: integer
doc: Back-end used in TREXIO. 0: HDF5, 1:Text
interface: ezfio, ocaml, provider
default: 0
[trexio_file]
type: character*(256)
doc: Name of the exported TREXIO file
interface: ezfio, ocaml, provider
default: None
[export_rdm]
type: logical
doc: If True, export two-body reduced density matrix
interface: ezfio, ocaml, provider
default: False
[export_ao_one_e_ints]
type: logical
doc: If True, export one-electron integrals in AO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_one_e_ints]
type: logical
doc: If True, export one-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False
[export_ao_two_e_ints]
type: logical
doc: If True, export two-electron integrals in AO basis
interface: ezfio, ocaml, provider
default: False
[export_ao_two_e_ints_cholesky]
type: logical
doc: If True, export Cholesky-decomposed two-electron integrals in AO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_two_e_ints]
type: logical
doc: If True, export two-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_two_e_ints_cholesky]
type: logical
doc: If True, export Cholesky-decomposed two-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False

6
src/trexio/README.rst Normal file
View File

@ -0,0 +1,6 @@
======
trexio
======
Module for handling TREXIO files.
See https://github.com/trex-coe/trexio

View File

@ -0,0 +1,7 @@
program export_trexio_prog
implicit none
read_wf = .True.
SOFT_TOUCH read_wf
call export_trexio
end

View File

@ -0,0 +1,604 @@
subroutine export_trexio
use trexio
implicit none
BEGIN_DOC
! Exports the wave function in TREXIO format
END_DOC
integer(trexio_t) :: f ! TREXIO file handle
integer(trexio_exit_code) :: rc
double precision, allocatable :: factor(:)
print *, 'TREXIO file : '//trim(trexio_filename)
print *, ''
call system('cp '//trim(trexio_filename)//' '//trim(trexio_filename)//'.bak')
if (backend == 0) then
f = trexio_open(trexio_filename, 'u', TREXIO_HDF5, rc)
else if (backend == 1) then
f = trexio_open(trexio_filename, 'u', TREXIO_TEXT, rc)
endif
if (f == 0_8) then
print *, 'Unable to open TREXIO file for writing'
print *, 'rc = ', rc
stop -1
endif
call ezfio_set_trexio_trexio_file(trexio_filename)
! ------------------------------------------------------------------------------
! Electrons
! ---------
print *, 'Electrons'
rc = trexio_write_electron_up_num(f, elec_alpha_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_electron_dn_num(f, elec_beta_num)
call trexio_assert(rc, TREXIO_SUCCESS)
! Nuclei
! ------
print *, 'Nuclei'
rc = trexio_write_nucleus_num(f, nucl_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_nucleus_charge(f, nucl_charge)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_nucleus_coord(f, nucl_coord_transp)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_nucleus_label(f, nucl_label, 32)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_nucleus_repulsion(f, nuclear_repulsion)
call trexio_assert(rc, TREXIO_SUCCESS)
! Pseudo-potentials
! -----------------
if (do_pseudo) then
print *, 'ECP'
integer :: num
num = 0
do k=1,pseudo_klocmax
do i=1,nucl_num
if (pseudo_dz_k(i,k) /= 0.d0) then
num = num+1
end if
end do
end do
do l=0,pseudo_lmax
do k=1,pseudo_kmax
do i=1,nucl_num
if (pseudo_dz_kl(i,k,l) /= 0.d0) then
num = num+1
end if
end do
end do
end do
integer, allocatable :: ang_mom(:), nucleus_index(:), power(:), lmax(:)
double precision, allocatable :: exponent(:), coefficient(:)
allocate(ang_mom(num), nucleus_index(num), exponent(num), coefficient(num), power(num), &
lmax(nucl_num) )
do i=1,nucl_num
lmax(i) = -1
do l=0,pseudo_lmax
do k=1,pseudo_kmax
if (pseudo_dz_kl_transp(k,l,i) /= 0.d0) then
lmax(i) = max(lmax(i), l)
end if
end do
end do
end do
j = 0
do i=1,nucl_num
do k=1,pseudo_klocmax
if (pseudo_dz_k_transp(k,i) /= 0.d0) then
j = j+1
ang_mom(j) = lmax(i)+1
nucleus_index(j) = i
exponent(j) = pseudo_dz_k_transp(k,i)
coefficient(j) = pseudo_v_k_transp(k,i)
power(j) = pseudo_n_k_transp(k,i)
end if
end do
do l=0,lmax(i)
do k=1,pseudo_kmax
if (pseudo_dz_kl_transp(k,l,i) /= 0.d0) then
j = j+1
ang_mom(j) = l
nucleus_index(j) = i
exponent(j) = pseudo_dz_kl_transp(k,l,i)
coefficient(j) = pseudo_v_kl_transp(k,l,i)
power(j) = pseudo_n_kl_transp(k,l,i)
end if
end do
end do
end do
lmax(:) = lmax(:)+1
rc = trexio_write_ecp_max_ang_mom_plus_1(f, lmax)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_z_core(f, int(nucl_charge_remove))
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_num(f, num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_ang_mom(f, ang_mom)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_exponent(f, exponent)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_coefficient(f, coefficient)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ecp_power(f, power)
call trexio_assert(rc, TREXIO_SUCCESS)
endif
! Basis
! -----
print *, 'Basis'
rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_prim_num(f, prim_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_shell_num(f, shell_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
call trexio_assert(rc, TREXIO_SUCCESS)
allocate(factor(shell_num))
if (ao_normalized) then
factor(1:shell_num) = shell_normalization_factor(1:shell_num)
else
factor(1:shell_num) = 1.d0
endif
rc = trexio_write_basis_shell_factor(f, factor)
call trexio_assert(rc, TREXIO_SUCCESS)
deallocate(factor)
rc = trexio_write_basis_shell_index(f, shell_index)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_exponent(f, prim_expo)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_basis_coefficient(f, prim_coef)
call trexio_assert(rc, TREXIO_SUCCESS)
allocate(factor(prim_num))
if (primitives_normalized) then
factor(1:prim_num) = prim_normalization_factor(1:prim_num)
else
factor(1:prim_num) = 1.d0
endif
rc = trexio_write_basis_prim_factor(f, factor)
call trexio_assert(rc, TREXIO_SUCCESS)
deallocate(factor)
! Atomic orbitals
! ---------------
print *, 'AOs'
rc = trexio_write_ao_num(f, ao_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ao_cartesian(f, 1)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ao_shell(f, ao_shell)
call trexio_assert(rc, TREXIO_SUCCESS)
integer :: i, pow0(3), powA(3), j, k, l, nz
double precision :: normA, norm0, C_A(3), overlap_x, overlap_z, overlap_y, c
nz=100
C_A(1) = 0.d0
C_A(2) = 0.d0
C_A(3) = 0.d0
allocate(factor(ao_num))
if (ao_normalized) then
do i=1,ao_num
l = ao_first_of_shell(ao_shell(i))
factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
enddo
else
factor(:) = 1.d0
endif
rc = trexio_write_ao_normalization(f, factor)
call trexio_assert(rc, TREXIO_SUCCESS)
deallocate(factor)
! One-e AO integrals
! ------------------
if (export_ao_one_e_ints) then
print *, 'AO one-e integrals'
rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ao_1e_int_kinetic(f,ao_kinetic_integrals)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_ao_1e_int_potential_n_e(f,ao_integrals_n_e)
call trexio_assert(rc, TREXIO_SUCCESS)
if (do_pseudo) then
rc = trexio_write_ao_1e_int_ecp(f, ao_pseudo_integrals_local + ao_pseudo_integrals_non_local)
call trexio_assert(rc, TREXIO_SUCCESS)
endif
rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
! Two-e AO integrals
! ------------------
if (export_ao_two_e_ints) then
print *, 'AO two-e integrals'
PROVIDE ao_two_e_integrals_in_map
integer(8), parameter :: BUFSIZE=100000_8
double precision :: eri_buffer(BUFSIZE), integral
integer(4) :: eri_index(4,BUFSIZE)
integer(8) :: icount, offset
double precision, external :: get_ao_two_e_integral
icount = 0_8
offset = 0_8
do l=1,ao_num
do k=1,ao_num
do j=l,ao_num
do i=k,ao_num
if (i==j .and. k<l) cycle
if (i<j) cycle
integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
if (integral == 0.d0) cycle
icount += 1_8
eri_buffer(icount) = integral
eri_index(1,icount) = i
eri_index(2,icount) = j
eri_index(3,icount) = k
eri_index(4,icount) = l
if (icount == BUFSIZE) then
rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
offset += icount
icount = 0_8
end if
end do
end do
end do
end do
if (icount >= 0_8) then
rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
end if
! Two-e AO integrals - Cholesky
! -----------------------------
integer(4) :: chol_index(3,BUFSIZE)
double precision :: chol_buffer(BUFSIZE)
if (export_ao_two_e_ints_cholesky) then
print *, 'AO two-e integrals Cholesky'
rc = trexio_write_ao_2e_int_eri_cholesky_num(f, cholesky_ao_num)
call trexio_assert(rc, TREXIO_SUCCESS)
icount = 0_8
offset = 0_8
do k=1,cholesky_ao_num
do j=1,ao_num
do i=1,ao_num
integral = cholesky_ao(i,j,k)
if (integral == 0.d0) cycle
icount += 1_8
chol_buffer(icount) = integral
chol_index(1,icount) = i
chol_index(2,icount) = j
chol_index(3,icount) = k
if (icount == BUFSIZE) then
rc = trexio_write_ao_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
offset += icount
icount = 0_8
end if
end do
end do
end do
if (icount > 0_8) then
rc = trexio_write_ao_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
end if
! Molecular orbitals
! ------------------
print *, 'MOs'
rc = trexio_write_mo_type(f, mo_label, len(trim(mo_label)))
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_mo_num(f, mo_num)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_mo_coefficient(f, mo_coef)
call trexio_assert(rc, TREXIO_SUCCESS)
if ( (trim(mo_label) == 'Canonical').and. &
(export_mo_two_e_ints_cholesky.or.export_mo_two_e_ints) ) then
rc = trexio_write_mo_energy(f, fock_matrix_diag_mo)
call trexio_assert(rc, TREXIO_SUCCESS)
endif
rc = trexio_write_mo_class(f, mo_class, len(mo_class(1)))
call trexio_assert(rc, TREXIO_SUCCESS)
! One-e MO integrals
! ------------------
if (export_mo_one_e_ints) then
print *, 'MO one-e integrals'
rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_mo_1e_int_potential_n_e(f,mo_integrals_n_e)
call trexio_assert(rc, TREXIO_SUCCESS)
if (do_pseudo) then
rc = trexio_write_mo_1e_int_ecp(f,mo_pseudo_integrals_local)
call trexio_assert(rc, TREXIO_SUCCESS)
endif
rc = trexio_write_mo_1e_int_core_hamiltonian(f,mo_one_e_integrals)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
! Two-e MO integrals
! ------------------
if (export_mo_two_e_ints) then
print *, 'MO two-e integrals'
PROVIDE mo_two_e_integrals_in_map
double precision, external :: mo_two_e_integral
icount = 0_8
offset = 0_8
do l=1,mo_num
do k=1,mo_num
do j=l,mo_num
do i=k,mo_num
if (i==j .and. k<l) cycle
if (i<j) cycle
integral = mo_two_e_integral(i,j,k,l)
if (integral == 0.d0) cycle
icount += 1_8
eri_buffer(icount) = integral
eri_index(1,icount) = i
eri_index(2,icount) = j
eri_index(3,icount) = k
eri_index(4,icount) = l
if (icount == BUFSIZE) then
rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
offset += icount
icount = 0_8
end if
end do
end do
end do
end do
if (icount > 0_8) then
rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
end if
! Two-e MO integrals - Cholesky
! -----------------------------
if (export_mo_two_e_ints_cholesky) then
print *, 'MO two-e integrals Cholesky'
rc = trexio_write_mo_2e_int_eri_cholesky_num(f, cholesky_ao_num)
call trexio_assert(rc, TREXIO_SUCCESS)
icount = 0_8
offset = 0_8
do k=1,cholesky_ao_num
do j=1,mo_num
do i=1,mo_num
integral = cholesky_mo(i,j,k)
if (integral == 0.d0) cycle
icount += 1_8
chol_buffer(icount) = integral
chol_index(1,icount) = i
chol_index(2,icount) = j
chol_index(3,icount) = k
if (icount == BUFSIZE) then
rc = trexio_write_mo_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
offset += icount
icount = 0_8
end if
end do
end do
end do
if (icount > 0_8) then
rc = trexio_write_mo_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
end if
! One-e RDM
! ---------
rc = trexio_write_rdm_1e(f,one_e_dm_mo)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_rdm_1e_up(f,one_e_dm_mo_alpha_average)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_rdm_1e_dn(f,one_e_dm_mo_beta_average)
call trexio_assert(rc, TREXIO_SUCCESS)
! Two-e RDM
! ---------
if (export_rdm) then
PROVIDE two_e_dm_mo
print *, 'Two-e RDM'
icount = 0_8
offset = 0_8
do l=1,mo_num
do k=1,mo_num
do j=1,mo_num
do i=1,mo_num
integral = two_e_dm_mo(i,j,k,l)
if (integral == 0.d0) cycle
icount += 1_8
eri_buffer(icount) = integral
eri_index(1,icount) = i
eri_index(2,icount) = j
eri_index(3,icount) = k
eri_index(4,icount) = l
if (icount == BUFSIZE) then
rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
offset += icount
icount = 0_8
end if
end do
end do
end do
end do
if (icount >= 0_8) then
rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
end if
end if
! ------------------------------------------------------------------------------
! Determinants
! ------------
integer*8, allocatable :: det_buffer(:,:,:)
double precision, allocatable :: coef_buffer(:,:)
integer :: nint
! rc = trexio_read_determinant_int64_num(f, nint)
! call trexio_assert(rc, TREXIO_SUCCESS)
nint = N_int
if (nint /= N_int) then
stop 'Problem with N_int'
endif
allocate ( det_buffer(nint, 2, BUFSIZE), coef_buffer(BUFSIZE, n_states) )
icount = 0_8
offset = 0_8
rc = trexio_write_state_num(f, n_states)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_set_state (f, 0)
call trexio_assert(rc, TREXIO_SUCCESS)
do k=1,n_det
icount += 1_8
det_buffer(1:nint, 1:2, icount) = psi_det(1:N_int, 1:2, k)
coef_buffer(icount,1:N_states) = psi_coef(k,1:N_states)
if (icount == BUFSIZE) then
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_determinant_list(f, offset, icount, det_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
do i=1,N_states
rc = trexio_set_state (f, i-1)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_determinant_coefficient(f, offset, icount, coef_buffer(1,i))
end do
rc = trexio_set_state (f, 0)
offset += icount
icount = 0_8
end if
end do
if (icount >= 0_8) then
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_determinant_list(f, offset, icount, det_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
do i=1,N_states
rc = trexio_set_state (f, i-1)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_write_determinant_coefficient(f, offset, icount, coef_buffer(1,i))
end do
rc = trexio_set_state (f, 0)
end if
deallocate ( det_buffer, coef_buffer )
rc = trexio_close(f)
call trexio_assert(rc, TREXIO_SUCCESS)
end
! -*- mode: f90 -*-

View File

@ -0,0 +1,79 @@
program import_determinants_ao
call run
end
subroutine run
use trexio
use map_module
implicit none
BEGIN_DOC
! Program to import determinants from TREXIO
END_DOC
integer(trexio_t) :: f ! TREXIO file handle
integer(trexio_exit_code) :: rc
integer :: m
double precision, allocatable :: coef_buffer(:,:)
integer*8 , allocatable :: det_buffer(:,:,:)
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
! Determinants
! ------------
integer :: nint, nstates
integer :: bufsize
rc = trexio_read_state_num(f, nstates)
call trexio_assert(rc, TREXIO_SUCCESS)
! rc = trexio_read_determinant_int64_num(f, nint)
! call trexio_assert(rc, TREXIO_SUCCESS)
nint = N_int
if (nint /= N_int) then
stop 'Problem with N_int'
endif
integer*8 :: offset, icount
rc = trexio_read_determinant_num(f, bufsize)
call trexio_assert(rc, TREXIO_SUCCESS)
print *, 'N_det = ', bufsize
allocate ( det_buffer(nint, 2, bufsize), coef_buffer(bufsize, n_states) )
offset = 0_8
icount = bufsize
rc = trexio_read_determinant_list(f, offset, icount, det_buffer)
call trexio_assert(rc, TREXIO_SUCCESS)
if (icount /= bufsize) then
print *, 'error: bufsize /= N_det: ', bufsize, icount
stop -1
endif
do m=1,nstates
rc = trexio_set_state(f, m-1)
call trexio_assert(rc, TREXIO_SUCCESS)
rc = trexio_read_determinant_coefficient(f, offset, icount, coef_buffer(1,m))
call trexio_assert(rc, TREXIO_SUCCESS)
if (icount /= bufsize) then
print *, 'error: bufsize /= N_det for state', m, ':', icount, bufsize
stop -1
endif
enddo
call save_wavefunction_general(bufsize,nstates,det_buffer,size(coef_buffer,1),coef_buffer)
end

View File

@ -0,0 +1,146 @@
program import_integrals_ao
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(f)
use trexio
use map_module
implicit none
BEGIN_DOC
! Program to import integrals from TREXIO
END_DOC
integer(trexio_t), intent(in) :: f ! TREXIO file handle
integer(trexio_exit_code) :: rc
integer ::i,j,k,l
integer(8) :: m, n_integrals
double precision :: integral
integer(key_kind), allocatable :: buffer_i(:)
real(integral_kind), allocatable :: buffer_values(:)
double precision, allocatable :: A(:,:)
double precision, allocatable :: V(:)
integer , allocatable :: Vi(:,:)
double precision :: s
if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
rc = trexio_read_nucleus_repulsion(f, s)
call trexio_assert(rc, TREXIO_SUCCESS)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here, rc
print *, 'Error reading nuclear repulsion'
stop -1
endif
call ezfio_set_nuclei_nuclear_repulsion(s)
call ezfio_set_nuclei_io_nuclear_repulsion('Read')
endif
! AO integrals
! ------------
allocate(A(ao_num, ao_num))
if (trexio_has_ao_1e_int_overlap(f) == TREXIO_SUCCESS) then
rc = trexio_read_ao_1e_int_overlap(f, A)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO overlap'
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A)
call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap('Read')
endif
if (trexio_has_ao_1e_int_kinetic(f) == TREXIO_SUCCESS) then
rc = trexio_read_ao_1e_int_kinetic(f, A)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO kinetic integrals'
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A)
call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
endif
! if (trexio_has_ao_1e_int_ecp(f) == TREXIO_SUCCESS) then
! rc = trexio_read_ao_1e_int_ecp(f, A)
! if (rc /= TREXIO_SUCCESS) then
! print *, irp_here
! print *, 'Error reading AO ECP local integrals'
! stop -1
! endif
! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A)
! call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo('Read')
! endif
if (trexio_has_ao_1e_int_potential_n_e(f) == TREXIO_SUCCESS) then
rc = trexio_read_ao_1e_int_potential_n_e(f, A)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO potential N-e integrals'
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A)
call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read')
endif
deallocate(A)
! AO 2e integrals
! ---------------
PROVIDE ao_integrals_map
integer*4 :: BUFSIZE
BUFSIZE=ao_num**2
allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
allocate(Vi(4,BUFSIZE), V(BUFSIZE))
integer*8 :: offset, icount
offset = 0_8
icount = BUFSIZE
rc = TREXIO_SUCCESS
do while (icount == size(V))
rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
do m=1,icount
i = Vi(1,m)
j = Vi(2,m)
k = Vi(3,m)
l = Vi(4,m)
integral = V(m)
call two_e_integrals_index(i, j, k, l, buffer_i(m) )
buffer_values(m) = integral
enddo
call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
offset = offset + icount
if (rc /= TREXIO_SUCCESS) then
exit
endif
end do
n_integrals = offset
call map_sort(ao_integrals_map)
call map_unique(ao_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
end

415
src/trexio/qp_import_trexio.py Executable file
View File

@ -0,0 +1,415 @@
#!/usr/bin/env python3
"""
convert TREXIO file to EZFIO
Usage:
qp_import_trexio [-o EZFIO_DIR] FILE
Options:
-o --output=EZFIO_DIR Produced directory
by default is FILE.ezfio
"""
import sys
import os
import trexio
import numpy as np
from functools import reduce
from ezfio import ezfio
from docopt import docopt
try:
QP_ROOT = os.environ["QP_ROOT"]
QP_EZFIO = os.environ["QP_EZFIO"]
except KeyError:
print("Error: QP_ROOT environment variable not found.")
sys.exit(1)
else:
sys.path = [QP_EZFIO + "/Python",
QP_ROOT + "/install/resultsFile",
QP_ROOT + "/install",
QP_ROOT + "/scripts"] + sys.path
def generate_xyz(l):
def create_z(x,y,z):
return (x, y, l-(x+y))
def create_y(accu,x,y,z):
if y == 0:
result = [create_z(x,y,z)] + accu
else:
result = create_y([create_z(x,y,z)] + accu , x, y-1, z)
return result
def create_x(accu,x,y,z):
if x == 0:
result = create_y([], x,y,z) + accu
else:
xnew = x-1
ynew = l-xnew
result = create_x(create_y([],x,y,z) + accu , xnew, ynew, z)
return result
result = create_x([], l, 0, 0)
result.reverse()
return result
def write_ezfio(trexio_filename, filename):
try:
trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_TEXT)
except:
trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
ezfio.set_file(filename)
ezfio.set_trexio_trexio_file(trexio_filename)
print("Nuclei\t\t...\t", end=' ')
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)
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)
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("Electrons\t...\t", end=' ')
try:
num_beta = trexio.read_electron_dn_num(trexio_file)
except:
num_beta = sum(charge)//2
try:
num_alpha = trexio.read_electron_up_num(trexio_file)
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)
print("OK")
print("Basis\t\t...\t", end=' ')
shell_num = 0
try:
basis_type = trexio.read_basis_type(trexio_file)
if basis_type.lower() not in ["gaussian", "slater"]:
raise TypeError
shell_num = trexio.read_basis_shell_num(trexio_file)
prim_num = trexio.read_basis_prim_num(trexio_file)
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
nucl_index = trexio.read_basis_nucleus_index(trexio_file)
exponent = trexio.read_basis_exponent(trexio_file)
coefficient = trexio.read_basis_coefficient(trexio_file)
shell_index = trexio.read_basis_shell_index(trexio_file)
ao_shell = trexio.read_ao_shell(trexio_file)
ezfio.set_basis_basis("Read from TREXIO")
ezfio.set_basis_shell_num(shell_num)
ezfio.set_basis_prim_num(prim_num)
ezfio.set_basis_shell_ang_mom(ang_mom)
ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
ezfio.set_basis_prim_expo(exponent)
ezfio.set_basis_prim_coef(coefficient)
nucl_shell_num = []
prev = None
m = 0
for i in ao_shell:
if i != prev:
m += 1
if prev is None or nucl_index[i] != nucl_index[prev]:
nucl_shell_num.append(m)
m = 0
prev = i
assert (len(nucl_shell_num) == nucl_num)
shell_prim_num = []
prev = shell_index[0]
count = 0
for i in shell_index:
if i != prev:
shell_prim_num.append(count)
count = 0
count += 1
prev = i
shell_prim_num.append(count)
assert (len(shell_prim_num) == shell_num)
ezfio.set_basis_shell_prim_num(shell_prim_num)
ezfio.set_basis_shell_index([x+1 for x in shell_index])
ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
shell_factor = trexio.read_basis_shell_factor(trexio_file)
prim_factor = trexio.read_basis_prim_factor(trexio_file)
print("OK")
except:
print("None")
ezfio.set_ao_basis_ao_cartesian(True)
print("AOS\t\t...\t", end=' ')
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)
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
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 = []
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)
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)]
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])
# 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")
# _
# |\/| _ _ |_) _. _ o _
# | | (_) _> |_) (_| _> | _>
#
print("MOS\t\t...\t", end=' ')
labels = { "Canonical" : "Canonical",
"RHF" : "Canonical",
"BOYS" : "Localized",
"ROHF" : "Canonical",
"UHF" : "Canonical",
"Natural": "Natural" }
try:
label = labels[trexio.read_mo_type(trexio_file)]
except:
label = "None"
ezfio.set_mo_basis_mo_label(label)
try:
clss = trexio.read_mo_class(trexio_file)
core = [ i for i in clss if i.lower() == "core" ]
inactive = [ i for i in clss if i.lower() == "inactive" ]
active = [ i for i in clss if i.lower() == "active" ]
virtual = [ i for i in clss if i.lower() == "virtual" ]
deleted = [ i for i in clss if i.lower() == "deleted" ]
except trexio.Error:
pass
try:
mo_num = trexio.read_mo_num(trexio_file)
ezfio.set_mo_basis_mo_num(mo_num)
MoMatrix = trexio.read_mo_coefficient(trexio_file)
ezfio.set_mo_basis_mo_coef(MoMatrix)
mo_occ = [ 0. for i in range(mo_num) ]
for i in range(num_alpha):
mo_occ[i] += 1.
for i in range(num_beta):
mo_occ[i] += 1.
ezfio.set_mo_basis_mo_occ(mo_occ)
except:
pass
print("OK")
print("Pseudos\t\t...\t", end=' ')
ezfio.set_pseudo_do_pseudo(False)
if trexio.has_ecp_ang_mom(trexio_file):
ezfio.set_pseudo_do_pseudo(True)
max_ang_mom_plus_1 = trexio.read_ecp_max_ang_mom_plus_1(trexio_file)
z_core = trexio.read_ecp_z_core(trexio_file)
ang_mom = trexio.read_ecp_ang_mom(trexio_file)
nucleus_index = trexio.read_ecp_nucleus_index(trexio_file)
exponent = trexio.read_ecp_exponent(trexio_file)
coefficient = trexio.read_ecp_coefficient(trexio_file)
power = trexio.read_ecp_power(trexio_file)
lmax = max( max_ang_mom_plus_1 ) - 1
ezfio.set_pseudo_pseudo_lmax(lmax)
ezfio.set_pseudo_nucl_charge_remove(z_core)
prev_center = None
ecp = {}
for i in range(len(ang_mom)):
center = nucleus_index[i]
if center != prev_center:
ecp[center] = { "lmax": max_ang_mom_plus_1[center],
"zcore": z_core[center],
"contr": {} }
for j in range(max_ang_mom_plus_1[center]+1):
ecp[center]["contr"][j] = []
ecp[center]["contr"][ang_mom[i]].append( (coefficient[i], power[i], exponent[i]) )
prev_center = center
ecp_loc = {}
ecp_nl = {}
kmax = 0
klocmax = 0
for center in ecp:
ecp_nl [center] = {}
for k in ecp[center]["contr"]:
if k == ecp[center]["lmax"]:
ecp_loc[center] = ecp[center]["contr"][k]
klocmax = max(len(ecp_loc[center]), klocmax)
else:
ecp_nl [center][k] = ecp[center]["contr"][k]
kmax = max(len(ecp_nl [center][k]), kmax)
ezfio.set_pseudo_pseudo_klocmax(klocmax)
ezfio.set_pseudo_pseudo_kmax(kmax)
pseudo_n_k = [[0 for _ in range(nucl_num)] for _ in range(klocmax)]
pseudo_v_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
pseudo_dz_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
pseudo_n_kl = [[[0 for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
pseudo_v_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
pseudo_dz_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
for center in ecp_loc:
for k in range( len(ecp_loc[center]) ):
v, n, dz = ecp_loc[center][k]
pseudo_n_k[k][center] = n
pseudo_v_k[k][center] = v
pseudo_dz_k[k][center] = dz
ezfio.set_pseudo_pseudo_n_k(pseudo_n_k)
ezfio.set_pseudo_pseudo_v_k(pseudo_v_k)
ezfio.set_pseudo_pseudo_dz_k(pseudo_dz_k)
for center in ecp_nl:
for l in range( len(ecp_nl[center]) ):
for k in range( len(ecp_nl[center][l]) ):
v, n, dz = ecp_nl[center][l][k]
pseudo_n_kl[l][k][center] = n
pseudo_v_kl[l][k][center] = v
pseudo_dz_kl[l][k][center] = dz
ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
print("OK")
def get_full_path(file_path):
file_path = os.path.expanduser(file_path)
file_path = os.path.expandvars(file_path)
return file_path
if __name__ == '__main__':
ARGUMENTS = docopt(__doc__)
FILE = get_full_path(ARGUMENTS['FILE'])
trexio_filename = FILE
if ARGUMENTS["--output"]:
EZFIO_FILE = get_full_path(ARGUMENTS["--output"])
else:
EZFIO_FILE = "{0}.ezfio".format(FILE)
write_ezfio(trexio_filename, EZFIO_FILE)
sys.stdout.flush()

View File

@ -0,0 +1,20 @@
BEGIN_PROVIDER [ character*(1024), trexio_filename ]
implicit none
BEGIN_DOC
! Name of the TREXIO file
END_DOC
character*(1024) :: prefix
trexio_filename = trexio_file
if (trexio_file == 'None') then
prefix = trim(ezfio_work_dir)//trim(ezfio_filename)
if (backend == 0) then
trexio_filename = trim(prefix)//'.h5'
else if (backend == 1) then
trexio_filename = trim(prefix)
endif
endif
END_PROVIDER

View File

@ -0,0 +1 @@
#include "trexio_f.f90"