Commit d840d505 authored by Sarai Dery Folkestad's avatar Sarai Dery Folkestad
Browse files

Merge branch 'speedup_cube_generation' into 'development' 

Speedup .cube generation by ~260x

See merge request !1543
parents 720df390 aa647f5b

changelog.md

+1 −0
Original line number Diff line number Diff line
@@ -243,6 +243,7 @@ some roots were converged and others below the linear dependence threshold. eT-p 
- Optimized convergence SC-QED-HF by second order in eta parameters. eT-program/eT!1498

- Optimized ``get_eri_1der_libcint`` to remove redundancy in derivative integrals. eT-program/eT!1511

- Optimized low memory CC2 Jacobian by replacing c1 transformations of the integrals, by c1 transformations of the Cholesky vectors. eT-program/eT!1518, eT-program/eT!1519, eT-program/eT!1520

- Optimized cube file generator. eT-program/eT!1543



### Documentation

- Updates to README. eT-program/eT!1247, eT-program/eT!1487

cmake/eT_sources.cmake

+1 −0
Original line number Diff line number Diff line
@@ -436,6 +436,7 @@ set(eT_fortran_sources 
   src/tools/density_storer_class.F90

   src/tools/string_utilities.F90

   src/tools/visualization_class.F90

   src/tools/sprintf_wrapper.c

   src/tools/z_matrix_tool_class.F90

   src/tools/electrostatics/point_charges_class.F90

   src/tools/electrostatics/point_charge_class.F90

src/memory/memory_manager_class.F90

+73 −2
Original line number Diff line number Diff line
@@ -137,13 +137,15 @@ module memory_manager_class 
                                  alloc_i32_3,     &

                                  alloc_l_1,       &

                                  alloc_l_2,       &

                                  alloc_i_2_ranges

                                  alloc_i_2_ranges,&

                                  alloc_char_1



      generic, public :: dealloc => dealloc_real,           &

                                    dealloc_complex,        &

                                    dealloc_64bit_integer,  &

                                    dealloc_32bit_integer,  &

                                    dealloc_logical

                                    dealloc_logical,        &

                                    dealloc_char





      procedure, private :: alloc_r_1
@@ -165,12 +167,14 @@ module memory_manager_class 
      procedure, private :: alloc_l_1

      procedure, private :: alloc_l_2

      procedure, private :: alloc_i_2_ranges

      procedure, private :: alloc_char_1



      procedure, private :: dealloc_real

      procedure, private :: dealloc_complex

      procedure, private :: dealloc_64bit_integer

      procedure, private :: dealloc_32bit_integer

      procedure, private :: dealloc_logical

      procedure, private :: dealloc_char



      procedure, private :: batch_setup_1

      procedure, private :: batch_setup_2
@@ -1239,6 +1243,38 @@ contains 
   end subroutine alloc_l_2





   subroutine alloc_char_1(mem, array, L, N)

      !!

      !! Written by Marcus T. Lexander, 2025

      !!

      !! Allocates an array of character(len=L)

      !!

      implicit none



      class(memory_manager) :: mem



      character(len=:), dimension(:), allocatable :: array



      integer, intent(in) :: L, N ! length of each string and number of strings



      integer(i64) :: size_array ! Total size of array (L*N)

      integer :: error = 0



      character(len=100) :: error_msg



      size_array = L*N



      allocate(character(len=L) :: array(N), stat=error, errmsg=error_msg)



      if (error .ne. 0) then

         call mem%print_allocation_error(size_array, error_msg)

      endif



      call mem%update_memory_after_alloc(size_array, 1)



   end subroutine alloc_char_1





   subroutine dealloc_real(mem, array)

      !!

      !! Written by Rolf H. Myhre and Alexander C. Paul, 2019-2023
@@ -1439,6 +1475,41 @@ contains 
   end subroutine dealloc_logical





   subroutine dealloc_char(mem, array)

      !!

      !! Written by Marcus T. Lexander, 2025

      !!

      !! Deallocates a character array and updates the available memory accordingly.

      !!

      implicit none



      class(memory_manager),                        intent(inout) :: mem

      character(len=:), dimension(..), allocatable, intent(inout) :: array



      integer(i64) :: size_array, char_len

      integer :: error = 0



      character(len=100) :: error_msg



      if(.not. allocated(array)) call output%error_msg("Trying to deallocate an unallocated array!")



      char_len = len(array)

      size_array = size(array)



      select rank(array)

         rank(1)

            deallocate(array, stat=error, errmsg=error_msg)

         rank default ! GCC 10.2 does not allow a call to output%error_msg here

            error stop 'Deallocation not implemented for character array of rank 2+'

      end select



      if (error .ne. 0) call mem%print_deallocation_error(size_array, error_msg)



      mem%available = mem%available + char_len*size_array



   end subroutine dealloc_char





   subroutine print_allocation_error(size_array, error_msg)

      !!

      !! Written by Alexander C. Paul, March 2020

src/tools/sprintf_wrapper.c

0 → 100644
+5 −0
Original line number Diff line number Diff line 
#include <stdio.h>



int sprintf_real_(char *buf, char *fmt, double *x) {

    return sprintf(buf, fmt, *x);

}

src/tools/visualization_class.F90

+42 −13
Original line number Diff line number Diff line
@@ -373,12 +373,14 @@ contains 
   subroutine write_vector_to_cube_visualization(plotter, ao, vector, file_name)

      !!

      !! Written by Tor S. Haugland, 2019

      !! Modified by Marcus T. Lexander, 2025

      !!

      !! Write vector of densities to a .cube file according to the specification,

      !!    http://paulbourke.net/dataformats/cube/

      !!

      use output_file_class, only: output_file

      use atomic_center_class, only: atomic_center

      use timings_class, only: timings



      implicit none
@@ -393,12 +395,22 @@ contains 
      real(dp) :: x_max, y_max, z_max

      real(dp) :: dx

      real(dp) :: volume, integral

      integer :: ix, iy, iz, i_center, n_printed

      logical :: newline

      integer :: ix, iy, iz, i_center, n_printed, print_buffer_len

      type(output_file) :: cube_file



      type(atomic_center) :: center



      character(len=:), dimension(:), allocatable :: printing_buffers

      integer, dimension(:), allocatable :: offsets



      character(len=100) :: timer_name

      type(timings) :: timer



      write(timer_name, '(a,a)') 'Writing .cube file: ', trim(file_name)

      timer = timings(timer_name, pl='normal')



      call timer%turn_on()



      ! .cube files are usually in Bohr

      x_min = plotter%x_min * angstrom_to_bohr

      y_min = plotter%y_min * angstrom_to_bohr
@@ -434,31 +446,48 @@ contains 
                           ints=[center%number_], reals=center%coordinates * angstrom_to_bohr)

      end do



      n_printed = 0

      newline = .false.

      ! One lane of the printing_buffers contains the printed text for all yz-coordinates

      ! given a single x coordinate.

      ! Currently each number is given 13 characters to be printed (12 for the number + 1 space or newline)

      print_buffer_len = plotter%n_y * plotter%n_z * 13



      call mem%alloc(printing_buffers, print_buffer_len, plotter%n_x)

      call mem%alloc(offsets, plotter%n_x)



      !$omp parallel do schedule(static) private(ix,iy,iz,n_printed)

      do ix=1, plotter%n_x

         offsets(ix) = 1

         n_printed = 0



         do iy=1, plotter%n_y

            do iz=1, plotter%n_z

               call cube_file%printf('m', "(e13.6)(x1)", fs='(a)', &

                                     reals=[vector(ix, iy, iz)], adv=.false.)

               newline = .false.

               call sprintf_real(printing_buffers(ix)(offsets(ix):), '%12.5e', vector(ix, iy, iz))

               offsets(ix) = offsets(ix) + 12



               n_printed = n_printed + 1

               if (n_printed == 6) then

                  call cube_file%printf('m', "", fs='(a)')

                  newline = .true.

                  printing_buffers(ix)(offsets(ix):offsets(ix)) = new_line(' ')

                  n_printed = 0

               else

                  printing_buffers(ix)(offsets(ix):offsets(ix)) = ' '

               end if

               offsets(ix) = offsets(ix) + 1

            enddo

         enddo

         if (.not. newline) then

            call cube_file%printf('m', "", fs='(a)')

            n_printed = 0

         end if

      enddo

      !$omp end parallel do



      do ix = 1, plotter%n_x

         write (cube_file%unit_, "(A)") printing_buffers(ix)(1:offsets(ix)-2)

      end do



      call mem%dealloc(printing_buffers)

      call mem%dealloc(offsets)



      call cube_file%close_()



      call timer%turn_off()



   end subroutine write_vector_to_cube_visualization