Add eshift to restart file to preserve continuity

.github/workflows/gfortran.yml

@@ -282,7 +282,7 @@ jobs:
       fail-fast: false
       matrix:
          gcc_v: [10]
-         mpich_v: ["3.4.3", "4.0.2"]
+         mpich_v: ["3.4.3", "4.0.3"]
     env:
       # To speed-up MPICH build
       CFLAGS: -O0

src/sh_integ.F90

@@ -454,21 +454,28 @@ end subroutine sh_set_energy_shift
    subroutine sh_write_wf(outunit)
       integer, intent(in) :: outunit
       integer :: ist1
+      write (outunit, '(A4, ES24.16E3)') 'E0= ', Eshift
       do ist1 = 1, nstate
          write (outunit, *) cel_re(ist1), cel_im(ist1)
       end do
    end subroutine sh_write_wf
 
    subroutine sh_read_wf(inunit)
       integer, intent(in) :: inunit
-      integer :: ist1
+      integer :: ist, iost
+      character(len=4) :: tmp
 
       if (.not. allocated(cel_re)) then
          allocate (cel_re(nstate), cel_im(nstate))
       end if
 
-      do ist1 = 1, nstate
-         read (inunit, *) cel_re(ist1), cel_im(ist1)
+      read (inunit, '(A4, ES24.16E3)', iostat=iost) tmp, Eshift
+      if (iost /= 0 .or. tmp /= 'E0= ') then
+         ! TODO: Write a more helpful error message.
+         call fatal_error(__FILE__, __LINE__, 'Invalid restart file.')
+      end if
+      do ist = 1, nstate
+         read (inunit, *) cel_re(ist), cel_im(ist)
       end do
    end subroutine sh_read_wf
 

src/surfacehop.F90

@@ -183,8 +183,8 @@ subroutine sh_init(x, y, z, vx, vy, vz)
       ! When restarting, initial SH WF was already read from the restart file
       if (irest == 0) then
          call sh_set_initialwf(istate)
+         call sh_set_energy_shift(en_array(1))
       end if
-      call sh_set_energy_shift(en_array(1))
       call sh_select_integrator(integ)
 
       entot0 = en_array(istate) + ekin_v(vx, vy, vz)

tests/SH_BAECK-AN/restart.xyz.ref

@@ -7,6 +7,7 @@
   1.1778599292239769E-003  0.0000000000000000E+000  0.0000000000000000E+000
  Coefficients for SH
            2
+E0=  2.2865474744001987E-003
    1.5068826045192000E-002  -1.0561205877374700E-003
   0.99459700118789496      -0.10270598969373472     
  Cumulative averages of various estimators

tests/SH_BUTCHER/hopgeom.3.2.6.xyz.ref

@@ -1,8 +1,8 @@
            6
 
- C   -3.6574528884969195E-003   2.4481075768171160E-002  -1.2398438630111645E-002
- N    1.6316824587872324E-002  -2.4633906933593016E-002   1.3252612775723083     
- H   0.80790415456768250       -3.8001865676746620E-002 -0.41109531238267050     
- H   0.79534332776060712        6.4856965131140865E-002   1.8506984897643315     
- H  -0.96480793921340524       0.12917535122661419        1.9223258122190650     
- H  -0.82470169291263140      -0.10576191436306318      -0.47268441499861041     
+C  -3.6574528884969195E-003  2.4481075768171160E-002 -1.2398438630111645E-002
+N   1.6316824587872324E-002 -2.4633906933593016E-002  1.3252612775723083E+000
+H   8.0790415456768250E-001 -3.8001865676746620E-002 -4.1109531238267050E-001
+H   7.9534332776060712E-001  6.4856965131140865E-002  1.8506984897643315E+000
+H  -9.6480793921340524E-001  1.2917535122661419E-001  1.9223258122190650E+000
+H  -8.2470169291263140E-001 -1.0576191436306318E-001 -4.7268441499861041E-001

tests/SH_BUTCHER/input.in.norestart

@@ -0,0 +1,39 @@
+&general
+irandom=   729375542,         ! random seed
+idebug=2
+pot='molpro'		! where do we obtain forces? options are:g09,orca,tera,turbo,molpro,orca,nab,harm,morse,guillot,2dho 
+mdtype='sh',		
+nwalk=1,		! number of random walkers 
+imini=0,		! equilibration period (should be at least 2000), or number of steps for minimization
+istage=0,		! staging transformation (1), without staging (0)
+nstep=20,		! number of steps
+dt=20.6697,			! timestep [au]
+irest=0,		! should we restart from restart.xyz? (ignoring mini.dat)
+
+nwrite=1,		! how often some output should be printed (estimators etc.)
+nwritex=1,		! how often should we print coordinates?
+nrest=1,		! how often we print restart files?
+nwritev=0,		! how often we print velocities?
+xyz_units='bohr'
+/
+
+
+&nhcopt
+temp=0.0,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
+inose=0,		! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2
+/
+
+&sh
+istate_init=3,		! initial electronic state
+nstate=3,		! number of electronic states
+substep=100,		! number of substeps for solving ESCH
+deltae=100.,		! maximum energy difference [eV], for which we calculate NA coupling
+integ='butcher',	! integrator for ESCHE:euler,butcher or rk4
+couplings='analytic',	! non-adiabatic coupling terms 'analytic', 'baeck-an', 'none'
+nohop=0,		! 1 - hopping not allowed (default=0)
+decoh_alpha=0.1		! parameter for decoherence correction
+popthr=0.01
+energydifthr=0.80
+energydriftthr=0.80
+correct_decoherence=.false.
+/