********> bugfix
cd $AMBERHOME
patch -p0 -N < bugfix

Author: Daniel R. Roe
Date: 2010-12-21
Programs: sander
Description: Fixes ref coords problem with hybrid REMD and ntr>0 

-----------------------------------------------------------------------------
diff --git src/sander/remd.f src/sander/remd.f
index c62a616..cbc4d27 100644
--- src/sander/remd.f
+++ src/sander/remd.f
@@ -147,6 +147,7 @@ integer, save :: incluster,nclust,nclustdih
 ! Hybrid REMD Temp. Storage
 _REAL_, dimension(:), allocatable, private :: hybrid_coord
 _REAL_, dimension(:), allocatable, private :: hybrid_force
+_REAL_, dimension(:), allocatable, private :: hybrid_refc
 
 ! H-remd temporary coordinates, force
 _REAL_, dimension(:), allocatable, private :: xtemp, ftemp
@@ -257,12 +258,13 @@ subroutine remd_setup(numexchg,hybridgb,numwatkeep,temp0,mxvar,natom)
       !  not known. amoeba may not work with hybrid remd.
       !   if (sanderrank==0) write(6,'(a)') "HYBRID REMD: Allocating memory."
       allocate( hybrid_coord(3*natom + mxvar), &
-                hybrid_force(3*natom + mxvar + 40), stat=ierror)
+                hybrid_force(3*natom + mxvar + 40), &
+                hybrid_refc(3*natom + mxvar), stat=ierror)
       REQUIRE( ierror == 0 )
       
       ! If an output file was specified for hybrid REMD stripped
       !  coords trajectory (-hybridtraj FILE), open it. Only masters
-      !  will writ to it.
+      !  will write to it.
       if (sanderrank==0) then
          if (remstripcoord(1:1).eq.' ') then
            write(6,'(a)') &
@@ -273,7 +275,7 @@ subroutine remd_setup(numexchg,hybridgb,numwatkeep,temp0,mxvar,natom)
             write(6,'(a,a)') &
                "HYBRID REMD: Opening hybrid stripped coord output: ",&
                remstripcoord
-            write(remlog_unit,'(a,a)') "#   hybdridtraj= ",remstripcoord
+            if (worldrank == 0) write(remlog_unit,'(a)') "# Writing stripped trajectories"
             call amopen(remstripcoord_unit,remstripcoord,'U','F','W')
             write(remstripcoord_unit,'(a80)') "Hybrid stripped coords"
          endif
@@ -369,6 +371,8 @@ subroutine remd_cleanup()
    REQUIRE(ier==0)
    if (allocated(hybrid_force)) deallocate(hybrid_force, stat=ier)
    REQUIRE(ier==0)
+   if (allocated(hybrid_refc)) deallocate(hybrid_refc, stat=ier)
+   REQUIRE(ier==0)
 
    if (allocated(saveene)) deallocate(saveene, stat=ier)
    REQUIRE(ier==0)
@@ -1167,7 +1171,7 @@ end subroutine load_reservoir_structure
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 subroutine hybrid_remd_ene( &
               x,ix,ih,ipairs,qsetup,                        &
-              numwatkeep,hybridgb,igb,nspm,t,temp0,                  &
+              numwatkeep,hybridgb,igb,ntr,nspm,t,temp0,                  &
               ntb,cut,                                               &
               ener,vir,do_list_update,nstep,nitp,nits,onefac,loutfm  ) 
    use state
@@ -1181,7 +1185,7 @@ subroutine hybrid_remd_ene( &
    character(len=4) ih(*)
    logical qsetup
 ! md.h
-   integer numwatkeep, hybridgb, igb, nspm
+   integer numwatkeep, hybridgb, igb, ntr, nspm
    _REAL_ t,temp0
 ! memory.h
 !   integer natom, nres
@@ -1200,7 +1204,7 @@ subroutine hybrid_remd_ene( &
 !  ntb, and cut during stripped coord call to force.
 ! DAN ROE: Is dynamic allocation the way to go?
 ! NOTE: This allocation currently won't work with amoeba.
-   integer hybrid_natom, hybrid_ntb, ier, i
+   integer hybrid_natom, ref_natom, hybrid_ntb, ier, i
    _REAL_ hybrid_cut
 
 !--------------------
@@ -1239,11 +1243,19 @@ subroutine hybrid_remd_ene( &
    do i=1, natom
       hybrid_coord(3*i - 2) = x(lcrd+3*i - 3)
       hybrid_coord(3*i - 1) = x(lcrd+3*i - 2)
-      hybrid_coord(3*i)     = x(lcrd+3*i - 1)
+      hybrid_coord(3*i    ) = x(lcrd+3*i - 1)
       hybrid_force(3*i - 2) = x(lforce+3*i - 3)
       hybrid_force(3*i - 1) = x(lforce+3*i - 2)
-      hybrid_force(3*i)     = x(lforce+3*i - 1)
+      hybrid_force(3*i    ) = x(lforce+3*i - 1)
    enddo
+   ! Store reference coordinates
+   if ( ntr == 1 ) then
+      do i=1, natom
+         hybrid_refc(3*i - 2) = x(lcrdr+3*i - 3)
+         hybrid_refc(3*i - 1) = x(lcrdr+3*i - 2)
+         hybrid_refc(3*i    ) = x(lcrdr+3*i - 1)  
+      enddo
+   endif
    hybrid_natom=natom
    hybrid_ntb=ntb
    hybrid_cut=cut
@@ -1253,8 +1265,16 @@ subroutine hybrid_remd_ene( &
    if (sanderrank==0) write(6,'(a)') "HYBRID REMD: Stripping waters"
    call stripwat(natom,hybrid_coord,ix(i02),ih(m02),x(lmass), &
                  ix(i70),nspm,nres,numwatkeep)
+   ! Strip reference coordinates. Strip actual array since it is
+   ! not passed into force. Will be restored from hybrid_refc after.
+   if ( ntr == 1 ) then
+      ref_natom = hybrid_natom
+      call stripwat(ref_natom,x(lcrdr),ix(i02),ih(m02),x(lmass), &
+                    ix(i70),nspm,nres,numwatkeep)
+   endif
    if (sanderrank==0) then
       write(6,'(a,i8)') "HYBRID REMD: New natom= ",natom
+      if (ntr == 1) write(6,'(a,i8)') "HYBRID REMD: New ref natom= ",ref_natom
       !write (6,*) "Post-strip coordinates: "
       !write (6,*) (xx(lcrd+i-1),i=1,10)
       ! DAN ROE: write to the stripped coordinate file
@@ -1262,7 +1282,7 @@ subroutine hybrid_remd_ene( &
       !                          "REMD ", repnum, mdloop, &
       !                          remstep, mytargettemp
       if (hybridwritetraj) &
-         call corpac(hybrid_coord,1,natom*3,remstripcoord_unit,loutfm)
+         call corpac(hybrid_coord,1,natom*3,remstripcoord_unit,.true.)
    endif
    call mpi_barrier(commworld, ier)
 
@@ -1286,7 +1306,7 @@ subroutine hybrid_remd_ene( &
 !            write (6,*) (xx(lforce+i-1),i=1,10)
       call prntmd(nstep,nitp,nits,t,ener,onefac,7,.false.)
    endif
-   if (sanderrank==0) write(6,'(a,f12.6,a,f6.2)') &
+   if (sanderrank==0) write(6,'(a,f13.4,a,f6.2)') &
       "HYBRID REMD: myEptot= ",ener%pot%tot," myTargetTemp= ",temp0
 
    ! 4- Restore original natom, igb, cut
@@ -1295,7 +1315,14 @@ subroutine hybrid_remd_ene( &
    igb=0
    ntb=hybrid_ntb
    cut=hybrid_cut
-
+   ! Restore reference coordinates
+   if ( ntr == 1 ) then
+      do i=1, natom
+         x(lcrdr+3*i - 3) = hybrid_refc(3*i - 2) 
+         x(lcrdr+3*i - 2) = hybrid_refc(3*i - 1) 
+         x(lcrdr+3*i - 1) = hybrid_refc(3*i    ) 
+      enddo
+   endif
 ! DAN ROE: Debug
 !         if (master) then
 !            write (6,*) "Restored coordinates: "
diff --git src/sander/runmd.f src/sander/runmd.f
index 5d35953..d1e191b 100644
--- src/sander/runmd.f
+++ src/sander/runmd.f
@@ -3312,7 +3312,7 @@ subroutine runmd(xx,ix,ih,ipairs,x,winv,amass,f, &
          ! This is a hybrid REMD run. Get energy of stripped system for next
          !  exchange.
          call hybrid_remd_ene(xx,ix,ih,ipairs,qsetup,      &
-                              numwatkeep,hybridgb,igb,nspm,t,temp0, &
+                              numwatkeep,hybridgb,igb,ntr,nspm,t,temp0, &
                               ntb,cut,                              &
                               ener,ener%vir,do_list_update,nstep,   &
                               nitp,nits,onefac,loutfm )
