--- AmberTools/src/sander/getcor_orig.F90	2021-01-31 16:05:08.834801880 +0000
+++ AmberTools/src/sander/getcor.F90	2021-02-01 21:22:32.247591490 +0000
@@ -5,7 +5,7 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+ Read initial coordinates (and velocities) from restart file
 #ifdef MPI
-subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0,writeflag,solvph,solve)
+subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0,writeflag,solvph,solve,stagid)
 #else
 subroutine getcor(nr,x,v,f,ntx,box,irest,tt,writeflag)
 #endif
@@ -32,11 +32,12 @@
    use sander_lib, only  : check_inpcrd_overflow, get_num_tokens
 #ifdef MPI
    use remd, only : rem
+   use sgld, only : isgld
 #endif
 
    implicit none
 #  include "ew_cntrl.h"
-   integer lun,nr,ntx,irest
+   integer lun,nr,ntx,irest,stagid
    _REAL_ x(*),v(*),f(*),box(3),tt
 #ifdef MPI
    _REAL_ temp0,solvph,solve
@@ -108,10 +109,19 @@
       ! in REM, overwrite temp0, solvph and/or solve if it's present in inpcrd/restrt
       ! DAN ROE: May want to just use temps in input file, so only 
       !          do this if a restart is requested.
-      if(rem == 1 .and. lremd_values(1) /= NO_INPUT_VALUE_FLOAT .and. irest==1) then
-         write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
-                          &netcdf inpcrd file'
-         temp0 = lremd_values(1)
+      if(rem == 1 .and. irest==1) then
+        if (isgld > 0) then
+          if (lremd_values(1) /= NO_INPUT_VALUE_FLOAT) then
+            write(6, '(a)') '| Correcting stagid using value from inpcrd file'
+            stagid = lremd_values(1)
+          else
+            stagid = 0
+          end if
+        else if (lremd_values(1) /= NO_INPUT_VALUE_FLOAT) then
+          write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
+            &netcdf inpcrd file'
+          temp0 = lremd_values(1)
+        end if
       else if(rem == 4 .and. lremd_values(1) /= NO_INPUT_VALUE_FLOAT .and. irest==1) then
          write(6, '(a)') '| Overwriting solvph from mdin with solvph from &
                           &netcdf inpcrd file'
@@ -123,10 +133,17 @@
       ! Multi-D REMD
       else if(rem == -1 .and. irest==1 .and. isRstValid) then
          do i = 1, remd_values_dim
-           if (lremd_types(i) == 1 .and. lremd_values(i) /= NO_INPUT_VALUE_FLOAT) then
-             write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
-                              &netcdf inpcrd file'
-             temp0 = lremd_values(i)
+           if(isgld > 0) then
+             if (lremd_types(i) == 1 .and. lremd_values(i) /= NO_INPUT_VALUE_FLOAT) then
+               write(6, '(a)') '| Correcting stagid using value from inpcrd file'
+               stagid = lremd_values(1)
+             else
+               stagid = 0
+             end if
+           else if (lremd_types(i) == 1 .and. lremd_values(i) /= NO_INPUT_VALUE_FLOAT) then
+               write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
+                                &netcdf inpcrd file'
+               temp0 = lremd_values(i)
            else if (lremd_types(i) == 4 .and. lremd_values(i) /= NO_INPUT_VALUE_FLOAT) then
              write(6, '(a)') '| Overwriting solvph from mdin with solvph from &
                               &netcdf inpcrd file'
@@ -171,19 +188,27 @@
       ! REMD simulations to set solvph to zero.
       call get_num_tokens(line_test, nwords)
       
-      9010 format(i5,2e15.7)
-      9011 format(i6,2e15.7)
-      9012 format(i7,2e15.7)
-      9013 format(i8,2e15.7)
-
+      9010 format(i5,2e15.7,2e15.7)
+      9011 format(i6,2e15.7,2e15.7)
+      9012 format(i7,2e15.7,2e15.7)
+      9013 format(i8,2e15.7,2e15.7)
 #ifdef MPI
       ! in REM, overwrite temp0 if it's present in inpcrd/restrt
       ! DAN ROE: May want to just use temps in input file, so only 
       !          do this if a restart is requested.
-      if(rem == 1 .and. lremd_values(1) > 0 .and. irest==1) then
-         write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
+      if(rem == 1 .and. irest==1) then
+        if(isgld > 0) then
+          if (lremd_values(1) > 0) then
+            write(6, '(a)') '| Correcting stagid using value from inpcrd file'
+            stagid = lremd_values(1)
+          else
+            stagid = 0
+          end if
+        else if (lremd_values(1) > 0) then
+          write(6, '(a)') '| Overwriting temp0 from mdin with temp0 from &
                           &netcdf inpcrd file'
-         temp0 = lremd_values(1)
+          temp0 = lremd_values(1)
+        end if
       else if(rem == 4 .and. nwords == 3 .and. irest==1) then
          write(6, '(a)') '| Overwriting solvph from mdin with solvph from &
                           &netcdf inpcrd file'
--- AmberTools/src/sander/remd_orig.F90	2020-04-28 02:19:46.000000000 +0100
+++ AmberTools/src/sander/remd.F90	2021-02-01 21:14:17.627615805 +0000
@@ -331,7 +331,7 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+ Set up REMD run - open log, set up temperature table, etc
 subroutine remd1d_setup(numexchg, hybridgb, numwatkeep, temp0, &
-                        mxvar, natom, ig, solvph, solve)
+                        mxvar, natom, ig, solvph, solve, stagid)
 
    use sgld, only : isgld, sgft, tempsg, sorttempsg, tempsglookup
 
@@ -346,6 +346,7 @@
    integer, intent(in) :: mxvar
    integer, intent(in) :: natom
    integer, intent(in) :: ig
+   integer, intent(in out) :: stagid
 
    _REAL_, intent(in) :: temp0
    _REAL_, intent(in) :: solvph
@@ -561,8 +562,13 @@
            ! Sort temperatures
            call sorttempsg(numreps, statetable, tsgtable, sgfttable)
            ! Determine this replca's ID
-           stagid = tempsglookup(numreps, temp0, tempsg, sgft, &
-                                 statetable, tsgtable, sgfttable)
+           if (stagid == 0) then
+             stagid = tempsglookup(numreps, temp0, &
+               tempsg, sgft, statetable, tsgtable, sgfttable)
+           else
+             tempsg = tsgtable(stagid)
+             sgft   = sgfttable(stagid)
+           end if
             ! Check to make sure the first replica has no guiding effect
            if(tsgtable(1)>0 .and. tsgtable(1) /= statetable(1) &
                             .and. sgfttable(1) /= 0) then
--- AmberTools/src/sander/sander_orig.F90	2021-01-31 17:20:16.864253844 +0000
+++ AmberTools/src/sander/sander.F90	2021-02-01 21:14:42.013663741 +0000
@@ -102,7 +102,8 @@
   ! Replica Exchange Molecular Dynamics
   use remd, only : rem, mdloop, remd1d_setup, remd_exchange, reservoir_remd_exchange, &
                    remd_cleanup, hremd_exchange, ph_remd_exchange, e_remd_exchange, &
-                   multid_remd_setup, multid_remd_exchange, setup_pv_correction, rremd
+                   multid_remd_setup, multid_remd_exchange, setup_pv_correction, &
+                   rremd, stagid
 #  ifdef VERBOSE_REMD
   use remd, only : repnum
 #  endif
@@ -431,7 +432,7 @@
           else
 #ifdef MPI
             call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, &
-                        box, irest, t, temp0, .FALSE., solvph, solve)
+                        box, irest, t, temp0, .FALSE., solvph, solve, stagid)
 #else
             call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, &
                         box, irest, t, .FALSE.)
@@ -632,7 +633,7 @@
 #ifdef LES
 #  ifdef MPI
         call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, box, irest, t, &
-                    temp0les, .TRUE., solvph, solve)
+                    temp0les, .TRUE., solvph, solve, stagid)
 #  else
         call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, box, irest, t, &
                     .TRUE.)
@@ -644,7 +645,7 @@
         else
 #  ifdef MPI
           call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, box, irest, &
-                      t, temp0, .TRUE., solvph, solve)
+                      t, temp0, .TRUE., solvph, solve, stagid)
 #  else
           call getcor(nr, x(lcrd), x(lvel), x(lforce), ntx, box, irest, t, &
                       .TRUE.)
@@ -1587,7 +1588,7 @@
 
             ! 1D REMD. Set up temptable, open remlog, etc.
             call remd1d_setup(numexchg, hybridgb, numwatkeep, &
-                              temp0, mxvar, natom, ig, solvph, solve)
+                        temp0, mxvar, natom, ig, solvph, solve, stagid)
           end if
           call setup_pv_correction(6, ntp, sanderrank)
 

