Store dimensions of orbital matrix and CI vectors in Checkpoint.txt (#54)

Previously, the Checkpoint.txt would dump the vectors only,
meaning that OpenFMS must know upon restart
what the lengths of these vectors are.

This is awkward since it requires "early" communication
with the electronic structure backend upon restart that
is not needed in a normal run, in particular, we must 
  1) ask backend for dimensions
  2) read checkpoint
  3) ask backend for electronic structure at restart geometry. 

Adding dimensions to Checkpoint.txt allows for a simpler sequence: 
  1) read checkpoint
  2) ask backend for electronic structure at restart geometry.

The code should be backwards compatible with pre-existing Checkpoint.txt files,
and doesn't modify the TeraChem initialization sequence.

Author: eirik-kjonstad

src/modules/RestartModule.f90

@@ -670,13 +670,30 @@ subroutine WriteElecStruc(ES, nf)
 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       type(T_ElecStruc), intent(in) :: ES
       integer(kind=DefInt), intent(in) :: nf
+      integer(kind=DefInt) :: rows, cols
 
 3     format(10(1x, es15.8))
 
-      write (nf, *) '# Orbitals'
-      write (nf, 3) ES%OldOrbitals
-      write (nf, *) '# CI vectors'
-      write (nf, 3) ES%OldCIVecs
+      if (allocated(ES%OldOrbitals)) then
+         rows = size(ES%OldOrbitals, dim=1)
+         cols = size(ES%OldOrbitals, dim=2)
+         write (nf, *) '# Orbitals', rows, cols
+         write (nf, 3) ES%OldOrbitals
+      else
+         write (nf, *) '# Orbitals', 0, 0
+         write (nf, *)
+      end if
+
+      if (allocated(ES%OldCIVecs)) then
+         rows = size(ES%OldCIVecs, dim=1)
+         cols = size(ES%OldCIVecs, dim=2)
+         write (nf, *) '# CI vectors', rows, cols
+         write (nf, 3) ES%OldCIVecs
+      else
+         write (nf, *) '# CI vectors', 0, 0
+         write (nf, *)
+      end if
+
       write (nf, *) '# Overlap matrix'
       write (nf, 3) ES%OverlapMatrix
       write (nf, *) '# TC Blob'
@@ -685,28 +702,176 @@ subroutine WriteElecStruc(ES, nf)
       write (nf, 3) ES%OldMSPT2C
       write (nf, *) '# Electronic Phases'
       write (nf, *) ES%ElecPhase
+
    end subroutine WriteElecStruc
 
    subroutine ReadElecStruc(ES, nf)
+
+      use ElecStrucModule, only: esNBasis, esLCIVec
       type(T_ElecStruc), intent(inout) :: ES
       integer(kind=DefInt), intent(in) :: nf
+      integer(kind=DefInt) :: ierr, rows, cols
+      character(len=500) :: header
 
 3     format(10(1x, es15.8))
 
-      read (nf, *)
-      read (nf, 3) ES%OldOrbitals
-      read (nf, *)
-      read (nf, 3) ES%OldCIVecs
-      read (nf, *)
+      read (nf, '(A)') header
+
+      ! Read orbitals
+      call parse_matrix_header(header, '# Orbitals', rows, cols)
+      call resolve_orbital_restart_dimensions(rows, cols)
+      call read_matrix_from_unit(ES%OldOrbitals, nf, rows, cols, 'orbital matrix')
+      call read_next_restart_header(nf, header)
+      if (rows > 0) esNBasis = rows
+
+      ! Read CI vectors
+      call parse_matrix_header(header, '# CI vectors', rows, cols)
+      call resolve_ci_restart_dimensions(rows, cols, size(ES%PotEn))
+      call read_matrix_from_unit(ES%OldCIVecs, nf, rows, cols, 'CI vector matrix')
+      call read_next_restart_header(nf, header)
+      if (cols > 0) esLCIVec = cols
+
+      call require_header(header, '# Overlap matrix')
       read (nf, 3) ES%OverlapMatrix
       read (nf, *)
       read (nf, 3) ES%OldBlob
       read (nf, *)
       read (nf, 3) ES%OldMSPT2C
       read (nf, *)
       read (nf, *) ES%ElecPhase
+
    end subroutine ReadElecStruc
 
+   subroutine read_next_restart_header(nf, header)
+      !!
+      !! Reads the next non-empty restart header
+      !!
+      integer(kind=DefInt), intent(in) :: nf
+      character(len=500), intent(out) :: header
+      integer(kind=DefInt) :: ierr
+
+      do
+         read (nf, '(A)') header
+         if (len_trim(header) > 0) exit
+      end do
+
+   end subroutine read_next_restart_header
+
+   subroutine resolve_orbital_restart_dimensions(rows, cols)
+      !!
+      !! Resolve old-orbital matrix dimensions from the restart header or from
+      !! an already-initialized backend, and check that both sources agree.
+      !!
+      use ElecStrucModule, only: esNBasis
+      integer(kind=DefInt), intent(inout) :: rows, cols
+
+      if (rows >= 0 .and. cols >= 0) then
+         if (rows == 0 .and. cols == 0) return
+         if (esNBasis > 0 .and. (rows /= esNBasis .or. cols /= esNBasis)) then
+            call FMS_DieError('Orbital matrix dimensions in restart file do not match initialized backend')
+         end if
+      else
+         if (esNBasis <= 0) then
+            call FMS_DieError('Restart orbital header must include dimensions when esNBasis is unknown')
+         end if
+         rows = esNBasis
+         cols = esNBasis
+      end if
+
+   end subroutine resolve_orbital_restart_dimensions
+
+   subroutine resolve_ci_restart_dimensions(rows, cols, num_states)
+      !!
+      !! Resolve CI vector dimensions from the restart header or from an
+      !! already-initialized backend, and check that both sources agree.
+      !!
+      use ElecStrucModule, only: esLCIVec
+      integer(kind=DefInt), intent(inout) :: rows, cols
+      integer(kind=DefInt), intent(in) :: num_states
+
+      if (rows >= 0 .and. cols >= 0) then
+         if (rows == 0 .and. cols == 0) return
+         if (rows /= num_states) then
+            call FMS_DieError('CI vector row count in restart file does not match number of states')
+         end if
+         if (esLCIVec > 0 .and. cols /= esLCIVec) then
+            call FMS_DieError('CI vector length in restart file does not match initialized backend')
+         end if
+      else
+         if (esLCIVec <= 0) then
+            call FMS_DieError('Restart CI vector header must include dimensions when esLCIVec is unknown')
+         end if
+         rows = num_states
+         cols = esLCIVec
+      end if
+
+   end subroutine resolve_ci_restart_dimensions
+
+   subroutine read_matrix_from_unit(matrix, nf, rows, cols, name)
+      !!
+      !! Reads named matrix from nf unit
+      !!
+      real(kind=DefReal), allocatable, intent(inout) :: matrix(:, :)
+      integer(kind=DefInt), intent(in) :: nf, rows, cols
+      character(len=*), intent(in) :: name
+      integer(kind=DefInt) :: ierr
+
+3     format(10(1x, es15.8))
+
+      if (rows < 0 .or. cols < 0) then
+         call FMS_DieError('Invalid dimensions for '//trim(name)//' in restart file')
+      end if
+
+      if (allocated(matrix)) deallocate (matrix)
+      if (rows > 0 .and. cols > 0) then
+         allocate (matrix(rows, cols))
+         read (nf, 3, iostat=ierr) matrix
+         if (ierr /= 0) call FMS_DieError('Error reading '//trim(name)//' from restart file')
+      end if
+
+   end subroutine read_matrix_from_unit
+
+   subroutine parse_matrix_header(header, label, rows, cols)
+      !!
+      !! Reads header with expected label (e.g., "# CI vectors") and extracts
+      !! the dimensions of the stored CI vectors according to the format
+      !! "# CI vectors [#rows] [#cols]".
+      !!
+      character(len=*), intent(in) :: header, label
+      integer(kind=DefInt), intent(out) :: rows, cols
+      integer(kind=DefInt) :: ios, start_pos
+      character(len=len(header)) :: adjusted_header
+
+      call require_header(header, label)
+
+      rows = -1
+      cols = -1
+      adjusted_header = adjustl(header)
+      start_pos = len_trim(label) + 1
+      if (len_trim(adjusted_header) >= start_pos) then
+         read (adjusted_header(start_pos:), *, iostat=ios) rows, cols
+         if (ios /= 0) then
+            rows = -1
+            cols = -1
+         end if
+      end if
+
+   end subroutine parse_matrix_header
+
+   subroutine require_header(header, label)
+      !!
+      !! Sanity check that expected header is actually == label
+      !!
+      character(len=*), intent(in) :: header, label
+
+      if (index(adjustl(header), trim(label)) /= 1) then
+
+         call FMS_DieError('Unsupported restart trajectory format: expected '''//trim(label)//'''')
+
+      end if
+
+   end subroutine require_header
+
 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 !           PARTICLE
 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

tests/AIMS-RESTART/Last_Bundle.txt.initial

@@ -57,9 +57,9 @@
   1.206920724811525E-002  0.000000000000000E+000  0.000000000000000E+000
   5.598543749720466E+00                       / Phase
        (0.98890207074517156,-1.64587765755770362E-003)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
@@ -107,9 +107,9 @@
   1.666180306845395E-002  0.000000000000000E+000  0.000000000000000E+000
   4.876812830655000E-01                       / Phase
        (-0.14543955145958318,3.02873312433818639E-002)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
@@ -157,9 +157,9 @@
   1.517517889199943E-002  0.000000000000000E+000  0.000000000000000E+000
   4.819255093038093E-01                       / Phase
        (-0.14498755228926799,3.00677562119372778E-002)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00

tests/AIMS-RESTART/Last_Bundle.txt.ref

@@ -57,9 +57,9 @@
   1.235177642280282E-002  0.000000000000000E+000  0.000000000000000E+000
   5.577528940734926E+00                       / Phase
        (0.98831589475674242,-1.79586289170714135E-003)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
@@ -107,9 +107,9 @@
   1.678084999429902E-002  0.000000000000000E+000  0.000000000000000E+000
   5.109671902184337E-01                       / Phase
        (-0.14899984083154516,3.20548667099393517E-002)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
@@ -157,9 +157,9 @@
   1.541008853267562E-002  0.000000000000000E+000  0.000000000000000E+000
   4.819255093038093E-01                       / Phase
        (-0.14498755228926799,3.00677562119372778E-002)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00

tests/ENERGY_CONS_REJECT/Last_Bundle.txt.ref

@@ -57,9 +57,9 @@
   2.497120962482854E-003  0.000000000000000E+000  0.000000000000000E+000
   6.271759009116956E+00                       / Phase
         (1.0000000000000000,-5.96311194867027439E-019)  / Amplitude
- # Orbitals
+ # Orbitals           0           0
 
- # CI vectors
+ # CI vectors           0           0
 
  # Overlap matrix
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00