Add TREXIO support

INSTALL.md

@@ -491,6 +491,14 @@ DeePMD-kit - Deep Potential Molecular Dynamics. Support for DeePMD-kit can be en
 - Add `-D__SMEAGOL` to DFLAGS, `-I$(LIBSMEAGOL_DIR)/obj` to FCFLAGS and
   `-L$(LIBSMEAGOL_DIR)/lib -lsmeagol` to LIBS
 
+### 2z. TREXIO (optional, unified computational chemistry format)
+
+TREXIO - Open-source file format and library. Support for TREXIO can be enabled via the flag
+`-D__TREXIO`.
+
+- TREXIO library can be downloaded from <https://github.com/trex-coe/trexio>
+- For more information see <https://trex-coe.github.io/trexio/index.html>.
+
 ## 3. Compile
 
 ### 3a. ARCH files
@@ -580,8 +588,9 @@ libraries (see 2.)
 - `-D__CRAY_PM_ACCEL_ENERGY` or `-D__CRAY_PM_ENERGY` switch on energy profiling on Cray systems
 - `-D__NO_ABORT` to avoid calling abort, but STOP instead (useful for coverage testing, and to avoid
   core dumps on some systems)
-- `-D__HDF5` enables hdf5 support. This is a hard dependency for SIRIUS, but can also be used by
-  itself to allow read/write functionalities of QCSchema files in the active space module.
+- `-D__HDF5` enables hdf5 support. This is a hard dependency for SIRIUS and TREXIO, but can also be
+  used by itself to allow read/write functionalities of QCSchema files in the active space module
+- `-D__TREXIO` enables TREXIO I/O support
 
 Features useful to deal with legacy systems
 

src/CMakeLists.txt

@@ -826,6 +826,8 @@ list(
   torch_api.F
   transport_env_types.F
   transport.F
+  trexio.F
+  trexio_utils.F
   uff_vdw_radii_table.F
   virial_methods.F
   voronoi_interface.F

src/aobasis/basis_set_types.F

@@ -614,12 +614,13 @@ END SUBROUTINE combine_basis_sets
 !> \param nshell_sum ...
 !> \param maxder ...
 !> \param short_kind_radius ...
+!> \param npgf_seg_sum number of primitives in "segmented contraction format"
 ! **************************************************************************************************
    SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, &
                                 nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, &
                                 m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, &
                                 last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, &
-                                npgf_sum, nshell_sum, maxder, short_kind_radius)
+                                npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
 
       ! Get informations about a Gaussian-type orbital (GTO) basis set.
 
@@ -646,6 +647,7 @@ SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radiu
                                                             nshell_sum
       INTEGER, INTENT(IN), OPTIONAL                      :: maxder
       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: short_kind_radius
+      INTEGER, INTENT(OUT), OPTIONAL                     :: npgf_seg_sum
 
       INTEGER                                            :: iset, nder
 
@@ -736,6 +738,12 @@ SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radiu
       END IF
       IF (PRESENT(npgf_sum)) npgf_sum = SUM(gto_basis_set%npgf)
       IF (PRESENT(nshell_sum)) nshell_sum = SUM(gto_basis_set%nshell)
+      IF (PRESENT(npgf_seg_sum)) THEN
+         npgf_seg_sum = 0
+         DO iset = 1, gto_basis_set%nset
+            npgf_seg_sum = npgf_seg_sum + gto_basis_set%npgf(iset)*gto_basis_set%nshell(iset)
+         END DO
+      END IF
 
    END SUBROUTINE get_gto_basis_set
 

src/cp2k_info.F

@@ -267,6 +267,10 @@ FUNCTION cp2k_flags() RESULT(flags)
       flags = TRIM(flags)//" hdf5"
 #endif
 
+#if defined(__TREXIO)
+      flags = TRIM(flags)//" trexio"
+#endif
+
 #if defined(__OFFLOAD_UNIFIED_MEMORY)
       flags = TRIM(flags)//" offload_unified_memory"
 #endif

src/input_cp2k_print_dft.F

@@ -617,6 +617,23 @@ SUBROUTINE create_print_dft_section(section)
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 
+      CALL section_create(subsection, __LOCATION__, name="TREXIO", &
+                          description="Write a TREXIO file to disk.", &
+                          n_keywords=1, n_subsections=0, repeats=.FALSE.)
+      CALL keyword_create(keyword, __LOCATION__, name="FILENAME", &
+                          description="Body of Filename for the trexio file.", &
+                          usage="FILENAME {name}", default_c_val="TREXIO", &
+                          type_of_var=char_t)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+      CALL keyword_create(keyword, __LOCATION__, name="CARTESIAN", &
+                          description="Store the MOs in the Cartesian basis instead of the default spherical basis.", &
+                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
       CALL section_create(subsection, __LOCATION__, name="GAPW", &
                           description="Controls the printing of some gapw related information (debug).", &
                           n_keywords=0, n_subsections=1, repeats=.FALSE.)

src/qs_kind_types.F

@@ -867,6 +867,7 @@ END SUBROUTINE get_qs_kind
 !> \param basis_rcut ...
 !> \param basis_type ...
 !> \param total_zeff_corr ... [SGh]
+!> \param npgf_seg total number of primitive GTOs in "segmented contraction format"
 ! **************************************************************************************************
    SUBROUTINE get_qs_kind_set(qs_kind_set, &
                               all_potential_present, tnadd_potential_present, gth_potential_present, &
@@ -875,7 +876,7 @@ SUBROUTINE get_qs_kind_set(qs_kind_set, &
                               ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, &
                               nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, &
                               basis_rcut, &
-                              basis_type, total_zeff_corr)
+                              basis_type, total_zeff_corr, npgf_seg)
 
       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
       LOGICAL, INTENT(OUT), OPTIONAL :: all_potential_present, tnadd_potential_present, &
@@ -889,6 +890,7 @@ SUBROUTINE get_qs_kind_set(qs_kind_set, &
       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: basis_rcut
       CHARACTER(len=*), OPTIONAL                         :: basis_type
       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: total_zeff_corr
+      INTEGER, INTENT(OUT), OPTIONAL                     :: npgf_seg
 
       CHARACTER(len=default_string_length)               :: my_basis_type
       INTEGER                                            :: ikind, imax, lmax_rho0_kind, &
@@ -944,6 +946,7 @@ SUBROUTINE get_qs_kind_set(qs_kind_set, &
          IF (PRESENT(lmax_rho0)) lmax_rho0 = 0
          IF (PRESENT(basis_rcut)) basis_rcut = 0.0_dp
          IF (PRESENT(total_zeff_corr)) total_zeff_corr = 0.0_dp
+         IF (PRESENT(npgf_seg)) npgf_seg = 0
 
          nkind = SIZE(qs_kind_set)
          DO ikind = 1, nkind
@@ -1223,6 +1226,13 @@ SUBROUTINE get_qs_kind_set(qs_kind_set, &
                lmax_rho0 = MAX(lmax_rho0, lmax_rho0_kind)
             END IF
 
+            IF (PRESENT(npgf_seg)) THEN
+               IF (ASSOCIATED(tmp_basis_set)) THEN
+                  CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_seg_sum=n)
+                  npgf_seg = npgf_seg + n*qs_kind_set(ikind)%natom
+               END IF
+            END IF
+
          END DO
       ELSE
          CPABORT("The pointer qs_kind_set is not associated")

src/qs_scf_post_gpw.F

@@ -194,6 +194,7 @@ MODULE qs_scf_post_gpw
    USE scf_control_types,               ONLY: scf_control_type
    USE stm_images,                      ONLY: th_stm_image
    USE transport,                       ONLY: qs_scf_post_transport
+   USE trexio_utils,                    ONLY: write_trexio
    USE virial_types,                    ONLY: virial_type
    USE voronoi_interface,               ONLY: entry_voronoi_or_bqb
    USE xray_diffraction,                ONLY: calculate_rhotot_elec_gspace,&
@@ -1685,7 +1686,7 @@ SUBROUTINE write_mo_dependent_results(qs_env, scf_env)
       CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_dependent_results'
 
       INTEGER                                            :: handle, homo, ispin, nmo, output_unit
-      LOGICAL                                            :: all_equal, do_kpoints
+      LOGICAL                                            :: all_equal, do_kpoints, explicit
       REAL(KIND=dp)                                      :: maxocc, s_square, s_square_ideal, &
                                                             total_abs_spin_dens
       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues, occupation_numbers
@@ -1712,7 +1713,8 @@ SUBROUTINE write_mo_dependent_results(qs_env, scf_env)
       TYPE(qs_rho_type), POINTER                         :: rho
       TYPE(qs_subsys_type), POINTER                      :: subsys
       TYPE(scf_control_type), POINTER                    :: scf_control
-      TYPE(section_vals_type), POINTER                   :: dft_section, input, sprint_section
+      TYPE(section_vals_type), POINTER                   :: dft_section, input, sprint_section, &
+                                                            trexio_section
 
       CALL timeset(routineN, handle)
 
@@ -1747,6 +1749,11 @@ SUBROUTINE write_mo_dependent_results(qs_env, scf_env)
       dft_section => section_vals_get_subs_vals(input, "DFT")
       IF (.NOT. qs_env%run_rtp) THEN
          CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
+         trexio_section => section_vals_get_subs_vals(dft_section, "PRINT%TREXIO")
+         CALL section_vals_get(trexio_section, explicit=explicit)
+         IF (explicit) THEN
+            CALL write_trexio(qs_env, trexio_section)
+         END IF
          IF (.NOT. do_kpoints) THEN
             CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv)
             CALL write_dm_binary_restart(mos, dft_section, ks_rmpv)

src/trexio.F

@@ -0,0 +1,17 @@
+!--------------------------------------------------------------------------------------------------!
+!   CP2K: A general program to perform molecular dynamics simulations                              !
+!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
+!                                                                                                  !
+!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
+!--------------------------------------------------------------------------------------------------!
+
+! **************************************************************************************************
+!> \brief Simple wrapper for the official TREXIO Fortran interface
+!> \par History
+!>      05.2024 created [SB]
+!> \author Stefano Battaglia
+! **************************************************************************************************
+
+#ifdef __TREXIO
+#include <trexio_f.f90>
+#endif