********>Bugfix 58:
Author: Bob Duke
Date: 12/08/2005

Programs: pmemd

Description: This patch is the equivalent of sander bugfixes 32, 54, and 56.

Fix: Apply the following patch to amber8/src/pmemd/src/runmd_cit.f90

------------------------------------------------------------------------------
*** runmd_cit.f90	2004-02-09 14:03:08.000000000 -0500
--- runmd_cit.f90    	2005-12-08 13:45:28.000000000 -0500
***************
*** 136,142 ****
    double precision      :: ecopy(51)
    double precision      :: vircopy(3)
  
-   equivalence (scaltp,     ener(5))
    equivalence (vol,        ener(10))
    equivalence (pres(1),    ener(11))
    equivalence (ekcmt(1),   ener(15))
--- 136,141 ----
***************
*** 468,473 ****
--- 467,493 ----
        end if
        call vrand_set_velocities(atm_cnt, vel, atm_mass_inv, temp0 * factt)
        if (belly) call bellyf_cit(atm_igroup, vel)
+ 
+       ! At this point in the code, the velocities lag the positions
+       ! by half a timestep.  If we intend for the velocities to be drawn
+       ! from a Maxwell distribution at the timepoint where the positions and
+       ! velocities are synchronized, we have to correct these newly
+       ! redrawn velocities by backing them up half a step using the
+       ! current force.
+       ! Note that this fix only works for Newtonian dynamics.
+ 
+       if (gammai .eq. 0.d0) then
+ #ifdef MPI
+         do atm_lst_idx = 1, my_atm_cnt
+           j = cit_my_atm_lst(atm_lst_idx)
+ #else
+         do j = 1, atm_cnt
+ #endif
+           wfac = atm_mass_inv(j) * half_dtx
+           vel(:, j) = vel(:, j) - frc(:, j) * wfac
+         end do
+       end if
+ 
      end if  ! (reset_velocities)
  
  ! Step 2: Do the velocity update:
  -------------------------------------------------------------------------------
  Temporary workarounds: none