From 49598822938da0ac9fbe9334f6f1d61d18de7f93 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 11 May 2023 22:48:48 +0200 Subject: [PATCH] Added TREXIO module --- src/trexio/EZFIO.cfg | 54 ++ src/trexio/README.rst | 6 + src/trexio/export_trexio.irp.f | 7 + src/trexio/export_trexio_routines.irp.f | 604 ++++++++++++++++++++ src/trexio/import_trexio_determinants.irp.f | 79 +++ src/trexio/import_trexio_integrals.irp.f | 146 +++++ src/trexio/qp_import_trexio.py | 415 ++++++++++++++ src/trexio/trexio_file.irp.f | 20 + src/trexio/trexio_module.F90 | 1 + 9 files changed, 1332 insertions(+) create mode 100644 src/trexio/EZFIO.cfg create mode 100644 src/trexio/README.rst create mode 100644 src/trexio/export_trexio.irp.f create mode 100644 src/trexio/export_trexio_routines.irp.f create mode 100644 src/trexio/import_trexio_determinants.irp.f create mode 100644 src/trexio/import_trexio_integrals.irp.f create mode 100755 src/trexio/qp_import_trexio.py create mode 100644 src/trexio/trexio_file.irp.f create mode 100644 src/trexio/trexio_module.F90 diff --git a/src/trexio/EZFIO.cfg b/src/trexio/EZFIO.cfg new file mode 100644 index 00000000..8606e908 --- /dev/null +++ b/src/trexio/EZFIO.cfg @@ -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 + diff --git a/src/trexio/README.rst b/src/trexio/README.rst new file mode 100644 index 00000000..7d7304c6 --- /dev/null +++ b/src/trexio/README.rst @@ -0,0 +1,6 @@ +====== +trexio +====== + +Module for handling TREXIO files. +See https://github.com/trex-coe/trexio diff --git a/src/trexio/export_trexio.irp.f b/src/trexio/export_trexio.irp.f new file mode 100644 index 00000000..3ae0dcb4 --- /dev/null +++ b/src/trexio/export_trexio.irp.f @@ -0,0 +1,7 @@ +program export_trexio_prog + implicit none + read_wf = .True. + SOFT_TOUCH read_wf + call export_trexio +end + diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f new file mode 100644 index 00000000..d69e7a70 --- /dev/null +++ b/src/trexio/export_trexio_routines.irp.f @@ -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= 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 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 -*- diff --git a/src/trexio/import_trexio_determinants.irp.f b/src/trexio/import_trexio_determinants.irp.f new file mode 100644 index 00000000..1759bb94 --- /dev/null +++ b/src/trexio/import_trexio_determinants.irp.f @@ -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 diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f new file mode 100644 index 00000000..9f9ad9d6 --- /dev/null +++ b/src/trexio/import_trexio_integrals.irp.f @@ -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 diff --git a/src/trexio/qp_import_trexio.py b/src/trexio/qp_import_trexio.py new file mode 100755 index 00000000..de8d1269 --- /dev/null +++ b/src/trexio/qp_import_trexio.py @@ -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() + diff --git a/src/trexio/trexio_file.irp.f b/src/trexio/trexio_file.irp.f new file mode 100644 index 00000000..c9897748 --- /dev/null +++ b/src/trexio/trexio_file.irp.f @@ -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 + + diff --git a/src/trexio/trexio_module.F90 b/src/trexio/trexio_module.F90 new file mode 100644 index 00000000..acd08492 --- /dev/null +++ b/src/trexio/trexio_module.F90 @@ -0,0 +1 @@ +#include "trexio_f.f90"