diff --git a/src/qmd/GNUmakefile b/src/qmd/GNUmakefile
index 32d95cadd0..bb2a03c556 100644
--- a/src/qmd/GNUmakefile
+++ b/src/qmd/GNUmakefile
@@ -4,7 +4,10 @@
 OBJ_OPTIMIZE =  qmd_init.o qmd_driver.o qmd_main.o qmd_input.o\
 		qmd_start.o task_qmd.o qmd_nucl_dt2.o qmd_nucl_dt.o\
 		qmd_gradient.o qmd_geom_store.o qmd_thermostat.o\
-		qmd_kinetic.o qmd_write_trj.o qmd_write_rst.o
+		qmd_kinetic.o qmd_write_trj.o qmd_write_rst.o\
+		qmd_ran1.o qmd_gasdev.o qmd_randomize.o\
+                qmd_nucl_pos.o qmd_nucl_vel.o qmd_read_rst.o\
+                qmd_gamdev.o qmd_sumno.o qmd_fix_momentum.o
 
 LIBRARY=libqmd.a
 
diff --git a/src/qmd/qmd_common.fh b/src/qmd/qmd_common.fh
index 5e0061e40c..ea6b80a839 100644
--- a/src/qmd/qmd_common.fh
+++ b/src/qmd/qmd_common.fh
@@ -7,12 +7,8 @@ c
       common/qmd_int/lgeom,igeom,jdebug
 c
 c     constants
-      double precision kb,au2fs
-      common/qmd_consts/kb,au2fs
-c
-c     berendsen thermostat
-      logical do_berendsen
-      double precision tau
+      double precision kb,au2fs,amu2au,au2ang
+      common/qmd_consts/kb,au2fs,amu2au,au2ang
 c
 c     nose-hoover thermostat parameters
       integer maxchain
@@ -22,9 +18,27 @@ c     nose-hoover thermostat parameters
       double precision v_nh(maxchain)    ! chain velocities
       double precision m_nh(maxchain)    ! chain masses
       double precision g_nh(maxchain)    ! bath couplings
+c
+c     berendsen thermostat
+      logical do_berendsen
+      double precision tau
+c
+c     langevin thermostat
+      logical do_langevin
+      double precision friction
+      integer idum
+c
+c     velocity rescaling
+      logical do_rescale
+c
+c     BDP rescaling
+      logical do_BDP
 c     
-      common/qmd_tstats/do_berendsen,tau,
-     &                      do_nosehoover,r_nh,v_nh,m_nh,g_nh
+      common/qmd_thermostats/do_berendsen,tau,
+     &                       do_nosehoover,r_nh,v_nh,m_nh,g_nh,
+     &                       do_langevin,friction,
+     &                       do_rescale,idum,
+     &                       do_BDP
 c
 c     temperature
       double precision targ_temp
diff --git a/src/qmd/qmd_driver.F b/src/qmd/qmd_driver.F
index 22313a2f9e..f56c67bc06 100644
--- a/src/qmd/qmd_driver.F
+++ b/src/qmd/qmd_driver.F
@@ -34,13 +34,18 @@ c
       integer nat                   ! number of atoms
       integer nactive               ! number of active atoms
       integer nstep_nucl           ! total nuclear steps
+      integer sstep_nucl           ! first nuclear step
       integer istep_nucl           ! nuclear md step
+      integer fstep                 ! how often to remove trans/rot
       double precision dt_nucl      ! nuclear md time step
+      double precision eltime       ! elapsed time (in fs)
       double precision esys         ! system energy (from task_gradient)
       double precision ekin         ! kinetic energy from the md part
       double precision etotal       ! total energy: esys + ekin
 c
       logical status
+      logical rst_exists
+      logical do_none
 c
       character*32 pname
       integer xyzfileunit
@@ -73,7 +78,7 @@ c     get rstfilename
       call util_file_name('qmdrst',.false.,.false.,rstfilename)
 c
 c     initialize parameters
-      call qmd_init(rtdb,nstep_nucl,dt_nucl)
+      call qmd_init(rtdb,nstep_nucl,dt_nucl,do_none,fstep)
 c
 c     create geometry object 
       if (.not. geom_create(geom,'geometry'))
@@ -106,7 +111,6 @@ c
       if(.not.ma_push_get(mt_dbl,3*nat,'grads',l_g,k_g))  ! grads
      $   call errquit(pname//'ma_push_get grads',3*nat,MA_ERR)
       call dfill(3*nat, 0.0d0, dbl_mb(k_g), 1)
-c
 c     get the active atoms (also need a check for linearity)
 c     need to implement constraints
       if (.not. ma_push_get(mt_log,nat,'active atoms',l_act,k_act)) ! active
@@ -123,25 +127,68 @@ c     get the masses, initial coordinates
      $   call errquit(pname//'geom_masses_get',geom,GEOM_ERR)
       if (.not.geom_cart_coords_get(geom,dbl_mb(k_r)))
      $   call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
+c
+      inquire(file=rstfilename,exist=rst_exists)
+      if (rst_exists) then
+        if (ga_nodeid().eq.0) then
+          write(luout,*) 'Restart file found'
+          write(luout,*) 'Reading positions, velocities, KE,'
+          write(luout,*) 'and last time step from file.'
+        end if
+c       read velocities from file
+        call qmd_read_rst(rstfilename,rstfileunit,nat,dbl_mb(k_v),
+     $         dbl_mb(k_r),ekin,sstep_nucl)
+c       remove global translations and rotations
+c       kinetic energy is the same on exit as was on entry
+        call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
+     $         ekin,.true.)
+c       update geom and store in rtdb
+        call qmd_geom_store(rtdb,geom,nat,dbl_mb(k_r))
+      else
+        if (ga_nodeid().eq.0) then
+          write(luout,*) 'No restart file found'
+          write(luout,*) 'Beginning with random velocities'
+        end if
+c       randomize velocities for the given temperature
+        call qmd_randomize(nat,dbl_mb(k_m),dbl_mb(k_v))
+        sstep_nucl=1
+c       remove global translations and rotations
+c       kinetic energy is NOT the same on exit as was on entry
+        call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
+     $         ekin,.false.)
+      end if
+c
+      call ga_sync()
+c
+      if (ga_nodeid().eq.0) then
+        final_temp = 2.0d0*ekin/kb/(1.d0*ndeg)
+        write(luout,120) "Current temp. (K):",final_temp
+      end if
 c
 c     perform energy gradient calculation on initial coordinates
       call qmd_gradient(rtdb,nat,dbl_mb(k_g),esys)
 c
-c     initialize velocity and acceleration 
-      call qmd_nucl_dt(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_g))
-c
-c     apply thermostat to initial velocities
-      call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
-c
 c     start nuclear molecular dynamics loop
-      do istep_nucl = 1,nstep_nucl
+      do istep_nucl = sstep_nucl,nstep_nucl
 c
-c       update position and velocity at half step: dt2
-        call qmd_nucl_dt2(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_r),
-     &      dbl_mb(k_v),dbl_mb(k_g))
+        if (mod(istep_nucl,fstep).eq.0) then
+c         remove global translations and rotations
+c         kinetic energy is the same on exit as was on entry
+          call qmd_fix_momentum(nat,dbl_mb(k_m),dbl_mb(k_v),dbl_mb(k_r),
+     $           ekin,.true.)
+        end if
+        if (.not.do_none) then
+c         apply thermostat to velocities
+          call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
+        end if
 c
-c       apply thermostat
-        call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
+c       update velocity at half step (dt/2)
+        call qmd_nucl_vel(nat,0.5d0*dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
+     &         dbl_mb(k_g))
+c
+c       update position at dt
+        call qmd_nucl_pos(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_r),
+     &        dbl_mb(k_v))
 c
 c       update geom and store in rtdb
         call qmd_geom_store(rtdb,geom,nat,dbl_mb(k_r))
@@ -149,17 +196,21 @@ c
 c       calculate energy and gradient
         call qmd_gradient(rtdb,nat,dbl_mb(k_g),esys)
 c
-c       update velocity at full step: dt
-        call qmd_nucl_dt(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
-     &             dbl_mb(k_g))
+c       update velocity at half step (dt/2)
+        call qmd_nucl_vel(nat,0.5d0*dt_nucl,dbl_mb(k_m),dbl_mb(k_v),
+     &         dbl_mb(k_g))
 c
-c       apply thermostat
-        call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
+        if (.not.do_none) then
+c         apply thermostat to velocities
+          call qmd_thermostat(nat,dt_nucl,dbl_mb(k_m),dbl_mb(k_v),ekin)
+        end if
 c
 c       persist updated velocity in geom
         if (.not.geom_vel_set(geom,dbl_mb(k_v)))
      $   call errquit(pname//'geom_vel_set',0,GEOM_ERR)
 c
+c       calculate current kinetic energy
+        call qmd_kinetic(nat,dbl_mb(k_m),dbl_mb(k_v),ekin)
 c       instantaneous total energy: instantaneous kinetic energy + system energy
         etotal = ekin + esys
 c
@@ -172,27 +223,28 @@ c
          write(luout,*)
          write(luout,*)
          call util_print_centered(LuOut,'QMD Run Information',20,.true.)
-         write(luout,*) "istep_nucl: ",istep_nucl
-         write(luout,*) "Time elapsed: ",istep_nucl*dt_nucl*au2fs
-         write(luout,*) "ekin: ",ekin
-         write(luout,*) "esys: ",esys
-         write(luout,*) "etotal: ",etotal
-         write(luout,*) "targ_temp: ",targ_temp
-         write(luout,*) "final_temp: ",final_temp
-         write(luout,*)
+         write(luout,110) "Current time step:",istep_nucl
+         eltime=dble(istep_nucl)*dt_nucl*au2fs
+         write(luout,130) "Time elapsed (fs):",eltime
+         write(luout,130) "Kinetic energy (a.u.):",ekin
+         write(luout,130) "Potential energy (a.u.):",esys
+         write(luout,130) "Total energy (a.u.):",etotal
+         write(luout,120) "Target temp. (K):",targ_temp
+         write(luout,120) "Current temp. (K):",final_temp
+
 c
 c        write trajectory
          call qmd_write_trj(xyzfilename,xyzfileunit,.true.,nat,
      &          istep_nucl,etotal,geom)
+c
+         call qmd_write_rst(rstfilename,rstfileunit,
+     &             nat,istep_nucl,ekin,geom)
         end if  ! ga_nodeid()
 c
       end do  ! istep_nucl
-c
-c     write restart information
-      if (ga_nodeid().eq.0) then
-         call qmd_write_rst(rstfilename,rstfileunit,
-     &             nat,istep_nucl-1,esys,ekin,geom)
-      end if
+  110 format(A25,I20)
+  120 format(A25,F20.2)
+  130 format(A25,F20.6)
 c
 c     clear memory
       if(.not.ma_pop_stack(l_act))! active
diff --git a/src/qmd/qmd_gradient.F b/src/qmd/qmd_gradient.F
index a989b1f217..b5b14a5c69 100644
--- a/src/qmd/qmd_gradient.F
+++ b/src/qmd/qmd_gradient.F
@@ -22,6 +22,8 @@ c
 c
       logical task_gradient
       external task_gradient
+      logical task_energy
+      external task_energy
 c
 c     preliminaries
       pname = "qmd_gradient: "
diff --git a/src/qmd/qmd_init.F b/src/qmd/qmd_init.F
index 747be0e92a..abea287c86 100644
--- a/src/qmd/qmd_init.F
+++ b/src/qmd/qmd_init.F
@@ -1,7 +1,7 @@
 c
 c     qmd_init: initialize parameters
 c     
-      subroutine qmd_init(rtdb,nstep_nucl,dt_nucl)
+      subroutine qmd_init(rtdb,nstep_nucl,dt_nucl,do_none,fstep)
 c
       implicit none
 c
@@ -22,6 +22,11 @@ c
       integer ao_bas_han            ! AO basis set handle
       integer nstep_nucl            ! nuclear steps
       integer nstep_elec            ! total electronic steps
+      integer states                ! number of electronic states
+      integer init_st               ! initially populated state
+      integer fstep                 ! how often to remove trans/rot
+      integer tvals(8)
+      integer svals
       double precision dt_nucl      ! nuclear md time step
       double precision dt_elec      ! electronic time step
       character*32 thermostat       ! thermostat name
@@ -29,6 +34,7 @@ c
       character*32 integrator_elec  ! electronic dynamics integrator
       integer ichain
       double precision nh_mass(2)   ! nose-hoover mass
+      logical do_none
 c
       character*32 pname
 c
@@ -45,57 +51,64 @@ c
        write(LuOut,*)
        write(LuOut,*)
        call util_print_centered
-     &  (LuOut,'QMD General Information',40,.true.)
+     &  (LuOut,'QMD Run Parameters',40,.true.)
        write(LuOut,*)
        write(LuOut,*)
       endif ! ga_nodeid
 c
 c     boltzmann constant in atomic units
-      kb=3.16679d-6
+      kb=3.1668114d-6
 c
 c     convert time in au to femtoseconds
-      au2fs=0.02419d0
+      au2fs=0.02418884326505d0
+c
+c     convert mass in amu to au
+      amu2au=1822.8884845d0
+c
+c     convert length in au to Angstrom
+      au2ang=0.52917721092d0
 c
 c     total nuclear steps
       if (.not.rtdb_get(rtdb,'qmd:nstep_nucl',mt_int,1,nstep_nucl))
      & call errquit(pname//'failed to read nstep_nucl',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "nstep_nucl: ",nstep_nucl
+      if (ga_nodeid().eq.0) 
+     &    write(luout,110) "No. of nuclear steps:",nstep_nucl
 c
 c     nuclear time step
       if (.not.rtdb_get(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
      & call errquit(pname//'failed to read dt_nucl',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "dt_nucl: ",dt_nucl
+      if (ga_nodeid().eq.0)
+     &    write(luout,120) "Nuclear time step:",dt_nucl
 c
 c     total electronic steps
       if(.not.rtdb_get(rtdb,'qmd:nstep_elec',mt_int,1,nstep_elec))
      & call errquit(pname//'failed to read nstep_elec',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "nstep_elec: ",nstep_elec
+c      if (ga_nodeid().eq.0) 
+c     &     write(luout,110) "No. of electronic steps:",nstep_elec
 c
 c     electronic time step
       if (.not.rtdb_get(rtdb,'qmd:dt_elec',mt_dbl,1,dt_elec))
      & call errquit(pname//'failed to read dt_elec',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "dt_elec: ",dt_elec
+c      if (ga_nodeid().eq.0)
+c     &     write(luout,*) "Electronic time step:",dt_elec
 c
 c     target temperature
       if (.not.rtdb_get(rtdb,'qmd:targ_temp',mt_dbl,1,targ_temp))
      & call errquit(pname//'failed to read targ_temp',0,RTDB_ERR)
       if (ga_nodeid().eq.0) 
-     &        write(luout,*) "target temperature: ",targ_temp
-c
-c     nose-hoover mass parameter
-      if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass))
-     & call errquit(pname//'failed to read nh_mass',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) 
-     &     write(luout,*) "nh_mass(1): ",nh_mass(1)
-      if (ga_nodeid().eq.0) 
-     &     write(luout,*) "nh_mass(2): ",nh_mass(2)
+     &        write(luout,120) "Target temp. (K):",targ_temp
 c
 c     thermostat
       do_berendsen = .false.
       do_nosehoover = .false.
+      do_langevin = .false.
+      do_BDP = .false.
+      do_rescale = .false.
+      do_none = .false.
       if (.not. rtdb_cget(rtdb,'qmd:thermostat',1,thermostat))
      $ call errquit(pname//'failed to read thermostat',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "thermostat: ",thermostat
+      if (ga_nodeid().eq.0)
+     $    write(luout,140) "Thermostat:",trim(thermostat)
       if (thermostat.eq.'berendsen') then
          do_berendsen = .true.
       else if (thermostat.eq.'nose-hoover') then
@@ -106,41 +119,86 @@ c     thermostat
             m_nh(ichain) = nh_mass(ichain) ! mass of chains
             g_nh(ichain) = 0.d0            ! bath coupling
          end do
+       else if (thermostat.eq.'langevin') then
+         do_langevin = .true.
+       else if (thermostat.eq.'rescale') then
+         do_rescale = .true.
+       else if (thermostat.eq.'BDP') then
+         do_BDP = .true.
+       else if (thermostat.eq.'none') then
+         do_none = .true.
       else
          if (ga_nodeid().eq.0) 
      &        write(luout,*) "unknown thermostat"
          if (ga_nodeid().eq.0) 
-     &        write(luout,*) "using default: berendsen thermostat"
-         do_berendsen = .true.
+     &        write(luout,*) "using default: NO thermostat"
+         do_none = .true.
       end if
 c
+c     nose-hoover mass parameter
+      if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass))
+     & call errquit(pname//'failed to read nh_mass',0,RTDB_ERR)
+      if (do_nosehoover) then
+       if (ga_nodeid().eq.0) 
+     &      write(luout,130) "NH mass(1):",nh_mass(1)
+       if (ga_nodeid().eq.0) 
+     &      write(luout,130) "NH mass(2):",nh_mass(2)
+      end if
+c
+c     langevin friction parameter
+      friction = 0.1d0
+      if (.not.rtdb_get(rtdb,'qmd:friction',mt_dbl,1,friction))
+     & friction = 0.1d0
+      if (do_langevin) then
+       if (ga_nodeid().eq.0) write(luout,130) "Friction:",friction
+      end if
+c
+c     berendsen tau parameter, also used for BDP
+      if (.not.rtdb_get(rtdb,'qmd:tau',mt_dbl,1,tau))
+     & call errquit(pname//'failed to read tau',0,RTDB_ERR)
+      if (do_berendsen.or.do_BDP) then
+       if (ga_nodeid().eq.0) write(luout,120) "Tau:",tau
+      end if
+c
+c     check dt_nucl against the berendsen tau
+      if (do_berendsen.and.(dt_nucl.ge.tau)) then
+        tau = 1.25d0*dt_nucl
+        if (ga_nodeid().eq.0) then
+         write(luout,*) "dt_nucl has to be smaller than tau"
+         write(luout,*) "adjusting tau to accomodate"
+         write(luout,*) "new tau: ", tau
+        end if
+      end if
+c
+c     random number seed
+      call date_and_time(values=tvals)
+c     current seconds x milliseconds + minutes
+      svals = tvals(7)*tvals(8)+tvals(6)
+      idum = -svals
+      if (.not.rtdb_get(rtdb,'qmd:rand_seed',mt_int,1,idum))
+     & idum = -svals ! default
+      if (ga_nodeid().eq.0) write(luout,110) "Random seed:",idum
+c
 c     nuclear md integrator
       if (.not. rtdb_cget(rtdb,'qmd:integrator_nucl',1,integrator_nucl))
      $ call errquit(pname//'failed to read integrator',0,RTDB_ERR)
       if (ga_nodeid().eq.0) 
-     & write(luout,*) "integrator_nucl: ",integrator_nucl
+     & write(luout,140) "Nuclear integrator:",trim(integrator_nucl)
       if (integrator_nucl.eq.'velocity-verlet') do_veloverlet = .true.
 c
-c     berendsen tau parameter
-      if (.not.rtdb_get(rtdb,'qmd:tau',mt_dbl,1,tau))
-     & call errquit(pname//'failed to read tau',0,RTDB_ERR)
-      if (ga_nodeid().eq.0) write(luout,*) "tau: ",tau
-c
-c     check dt_nucl against the berendsen tau
-      if (do_berendsen.and.(dt_nucl.ge.tau)) then
-        dt_nucl = 0.75d0*tau 
-        if (ga_nodeid().eq.0) then
-         write(luout,*) "dt_nucl has to be lesser than tau"
-         write(luout,*) "adjusting dt_nucl to accomodate"
-         write(luout,*) "new dt_nucl: ", dt_nucl
-        end if
-      end if
-c
 c     is linear molecule ?
       if (.not.rtdb_get(rtdb,'qmd:linear',mt_log,1,do_linear))
      &     do_linear = .false.
 c
+c     remove translations and rotations every fsteps
+      if (.not.rtdb_get(rtdb,'qmd:fstep',mt_int,1,fstep))
+     $     fstep = 1000
       if (ga_nodeid().eq.0) call util_flush(LuOut)
+c
+  110 format(A25,I20)
+  120 format(A25,F20.2)
+  130 format(A25,F20.6)
+  140 format(A25,A20)
 c
       return
       end
diff --git a/src/qmd/qmd_input.F b/src/qmd/qmd_input.F
index c4e812f520..9d1c28cfff 100644
--- a/src/qmd/qmd_input.F
+++ b/src/qmd/qmd_input.F
@@ -24,7 +24,10 @@ c
       integer nstep_nucl
       double precision dt_elec
       integer nstep_elec
+      integer fstep
       double precision targ_temp
+      double precision friction
+      integer idum
       double precision tau
       double precision nh_mass(2)
       character*32 test
@@ -42,7 +45,7 @@ c     Set default values
 c     ------------------
 c
 c     Nuclear MD time step: dt_nucl
-      dt_nucl = 1.d0
+      dt_nucl = 10.d0
       if (.not.rtdb_put(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
      & call errquit(pname//'failed to write dt_nucl',0,RTDB_ERR)
 c
@@ -67,12 +70,12 @@ c     Target nuclear MD temperature: targ_temp (Kelvin)
      & call errquit(pname//'failed to write targ_temp',0,RTDB_ERR)
 c
 c     Nuclear thermostat
-      thermostat='berendsen'
+      thermostat='none'
       if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
      & call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
 c
-c     Berendsen tau parameter
-      tau=0.1d0
+c     Berendsen tau parameter, also used for BDP
+      tau=1.d3
       if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
      & call errquit(pname//'failed to write tau',0,RTDB_ERR)
 c
@@ -93,6 +96,11 @@ c     Is linear molecule ?
       if (.not.rtdb_put(rtdb,'qmd:linear',mt_log,1,do_linear))
      & call errquit(pname//'failed to write do_linear',0,RTDB_ERR)
 c
+c     Remove translations/rotations every fsteps
+      fstep=1000
+      if (.not.rtdb_put(rtdb,'qmd:fstep',mt_int,1,fstep))
+     & call errquit(pname//'failed to write fstep',0,RTDB_ERR)
+c
 c     ----------
 c     Read input
 c     ----------
@@ -106,8 +114,8 @@ c     Nuclear MD time step: dt_nucl in femtoseconds
       if (inp_compare(.false.,test,'dt_nucl')) then
        if (.not.inp_f(dt_nucl)) then
           write(luout,*) 'dt_nucl: value not found; ',
-     1      'default value of 1.0 au will be used'
-          dt_nucl = 1.d0
+     1      'default value of 10.0 au will be used'
+          dt_nucl = 10.d0
        endif
        if (.not.rtdb_put(rtdb,'qmd:dt_nucl',mt_dbl,1,dt_nucl))
      &  call errquit(pname//'failed to write dt_nucl',0,RTDB_ERR)
@@ -169,10 +177,46 @@ c     Nuclear thermostat: Berendsen
           thermostat='berendsen'
        if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
      1  call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
-       if (.not.inp_f(tau)) tau = 0.1d0
+       if (.not.inp_f(tau)) tau = 1.d3
        if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
      &  call errquit(pname//'failed to write tau',0,RTDB_ERR)
 c
+c     Nuclear thermostat: BDP (Bussi, Donadio, Parrinello 2007)
+      else if (inp_compare(.false.,test,'BDP')) then
+          thermostat='BDP'
+       if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
+     1  call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
+       if (.not.inp_f(tau)) tau = 1.d2
+       if (.not.rtdb_put(rtdb,'qmd:tau',mt_dbl,1,tau))
+     &  call errquit(pname//'failed to write tau',0,RTDB_ERR)
+c
+c     Nuclear thermostat: Langevin (Bussi, Parrinello 2007))
+      else if (inp_compare(.false.,test,'langevin')) then
+          thermostat='langevin'
+       if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
+     1  call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
+c
+c     Nuclear thermostat: Velocity rescaling
+      else if (inp_compare(.false.,test,'rescale')) then
+          thermostat='rescale'
+       if (.not.rtdb_cput(rtdb,'qmd:thermostat',1,thermostat))
+     1  call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
+c
+c     Friction
+      else if (inp_compare(.false.,test,'friction')) then
+       if (.not.inp_f(friction)) friction = 0.1d0
+       if (.not.rtdb_put(rtdb,'qmd:friction',mt_dbl,1,friction))
+     &  call errquit(pname//'failed to write friction',0,
+     &     RTDB_ERR)
+c
+c     Random number seed
+      else if (inp_compare(.false.,test,'rand_seed')) then
+       if (.not.inp_i(idum)) idum = 3989
+       idum = sign(idum,-1)   ! seed needs to be negative
+       if (.not.rtdb_put(rtdb,'qmd:rand_seed',mt_int,1,idum))
+     &  call errquit(pname//'failed to write rand_seed',0,
+     &     RTDB_ERR)
+c
 c     Nuclear MD integrator
       else if (inp_compare(.false.,test,'velocity-verlet')) then
           integrator_nucl='velocity-verlet'
@@ -185,6 +229,12 @@ c     Is linear molecule ?
       if (.not.rtdb_put(rtdb,'qmd:linear',mt_log,1,do_linear))
      &  call errquit(pname//'failed to write do_linear',0,RTDB_ERR)
 c
+c     Remove translations and rotations every fsteps
+      else if (inp_compare(.false.,test,'fstep')) then
+       if (.not.inp_i(fstep)) fstep=1000
+       if (.not.rtdb_put(rtdb,'qmd:fstep',mt_int,1,fstep))
+     &  call errquit(pname//'failed to write fstep',0,RTDB_ERR)
+c
 c     Print
       else if (inp_compare(.false.,test,'print')) then
         call util_print_input(rtdb,'qmd')
diff --git a/src/qmd/qmd_kinetic.F b/src/qmd/qmd_kinetic.F
index 7d9892f093..9babb22ece 100644
--- a/src/qmd/qmd_kinetic.F
+++ b/src/qmd/qmd_kinetic.F
@@ -1,23 +1,25 @@
 c
 c     calculate kinetic energy 
 c
-      subroutine qmd_kinetic(nat,m,v,ekin)
+      subroutine qmd_kinetic(n,m,v,ekin)
 c
       implicit none
 c
-      integer nat
-      double precision m(nat)
-      double precision v(3,nat)
+#include "qmd_common.fh"
+c
+      integer n
+      double precision m(n)
+      double precision v(3,n)
       double precision ekin
 c
       integer i
 c
 c     total kinetic energy
       ekin = 0.d0
-      do i = 1,nat
+      do i = 1,n
        ekin = ekin +  m(i)*(v(1,i)*v(1,i)+v(2,i)*v(2,i)+v(3,i)*v(3,i))
       end do
-      ekin = 0.5d0*ekin
+      ekin = 0.5d0*ekin*amu2au
 c
       return
       end
diff --git a/src/qmd/qmd_nucl_dt.F b/src/qmd/qmd_nucl_dt.F
index 6cfbc2b938..f25c82bacd 100644
--- a/src/qmd/qmd_nucl_dt.F
+++ b/src/qmd/qmd_nucl_dt.F
@@ -1,36 +1,38 @@
 c
-c     qmd_nucl_dt: velocity at full dt
+c     Update velocity and positions
 c
-      subroutine qmd_nucl_dt(nat,dt,m,v,g)
+      subroutine qmd_nucl_dt(nat,dt,m,r,v,g)
 c
       implicit none
 c
       integer nat                ! number of atoms
       double precision dt        ! step size
       double precision m(nat)    ! mass
+      double precision r(3,nat)  ! coordinates
       double precision v(3,nat)  ! velocities
       double precision g(3,nat)  ! gradients
 c
+      integer i
       double precision a(3,nat)  ! acceleration
-      integer i,j
-      double precision dt2
 c
       character*30 pname
 c
 c     preliminaries
-      pname = "qmd_nucl_dt: "
-      dt2 = 0.5d0*dt
+      pname='qmd_nucl_dt: '
 c
-c     acceleration (a) and velocities (v) at dt
-      do i=1,nat
+c     calculate acceleration, update positions (r) and velocities (v) at dt/2
+      do i = 1,nat
         a(1,i) = -g(1,i)/m(i)
-        v(1,i) = v(1,i) + a(1,i)*dt2
+        v(1,i) = v(1,i) + a(1,i)*0.5d0*dt
+        r(1,i) = r(1,i) + v(1,i)*dt 
 c
         a(2,i) = -g(2,i)/m(i)
-        v(2,i) = v(2,i) + a(2,i)*dt2
+        v(2,i) = v(2,i) + a(2,i)*0.5d0*dt
+        r(2,i) = r(2,i) + v(2,i)*dt
 c
         a(3,i) = -g(3,i)/m(i)
-        v(3,i) = v(3,i) + a(3,i)*dt2
+        v(3,i) = v(3,i) + a(3,i)*0.5d0*dt
+        r(3,i) = r(3,i) + v(3,i)*dt
       end do
 c
       return
diff --git a/src/qmd/qmd_thermostat.F b/src/qmd/qmd_thermostat.F
index fcc56c9218..2a637be2e7 100644
--- a/src/qmd/qmd_thermostat.F
+++ b/src/qmd/qmd_thermostat.F
@@ -1,7 +1,7 @@
 c
 c     qmd_thermostat: apply thermostat
 c
-      subroutine qmd_thermostat(nat,dt,m,v,ekin)
+      subroutine qmd_thermostat(n,dt,m,v,ekin)
 c
       implicit none
 c
@@ -12,11 +12,12 @@ c
 c
 #include "qmd_common.fh"
 c
-      integer nat                ! number of atoms
-      double precision dt        ! step size
-      double precision m(nat)    ! mass
-      double precision v(3,nat)  ! velocities
-      double precision ekin      ! instantaneous kinetic energy
+      integer n                ! number of atoms
+      double precision dt      ! step size
+      double precision m(n)    ! mass
+      double precision v(3,n)  ! velocities
+      double precision ekin    ! instantaneous kinetic energy
+      double precision ekin2
 c
       character*30 pname
 c
@@ -24,13 +25,19 @@ c     preliminaries
       pname='qmd_thermostat: '
 c
 c     calculate kinetic energy
-      call qmd_kinetic(nat,m,v,ekin)
+      call qmd_kinetic(n,m,v,ekin)
 c
 c     apply thermostat
       if (do_nosehoover) then
-         call qmd_nosehoover(nat,dt,v,ekin)
+         call qmd_nosehoover(n,dt,m,v,ekin)
       else if (do_berendsen) then
-         call qmd_berendsen(nat,dt,v,ekin)
+         call qmd_berendsen(n,dt,m,v,ekin)
+      else if (do_langevin) then
+         call qmd_langevin(n,dt,m,v,ekin)
+      else if (do_rescale) then
+         call qmd_rescale(n,m,v,ekin,ekin2)
+      else if (do_BDP) then
+         call qmd_BDP(n,dt,m,v,ekin)
       else
          call errquit(pname//'unknown thermostat',0,RTDB_ERR)
       end if
@@ -40,7 +47,7 @@ c
 c
 c     qmd_nosehoover: nose-hoover thermostat
 c
-      subroutine qmd_nosehoover(nat,dt,v,ekin)
+      subroutine qmd_nosehoover(n,dt,m,v,ekin)
 c
       implicit none
 c
@@ -51,9 +58,10 @@ c
 c
 #include "qmd_common.fh"
 c
-      integer nat
+      integer n
       double precision dt
-      double precision v(3,nat)
+      double precision m(n)    ! mass
+      double precision v(3,n)
       double precision ekin
 c
       double precision dt2,dt4,dt8
@@ -102,7 +110,7 @@ c
       v_nh(2) = v_nh(2) + g_nh(2)*dt4
 c
 c     scale velocity
-      do i = 1,nat
+      do i = 1,n
         v(1,i) = scalefac*v(1,i)
         v(2,i) = scalefac*v(2,i)
         v(3,i) = scalefac*v(3,i)
@@ -113,7 +121,7 @@ c
 c
 c     qmd_berendsen: berendsen thermostat
 c
-      subroutine qmd_berendsen(nat,dt,v,ekin)
+      subroutine qmd_berendsen(n,dt,m,v,ekin)
 c
       implicit none
 c
@@ -124,9 +132,10 @@ c
 c
 #include "qmd_common.fh"
 c
-      integer nat
+      integer n
       double precision dt
-      double precision v(3,nat)
+      double precision m(n)    ! mass
+      double precision v(3,n)
       double precision ekin
 c
       integer i,j
@@ -142,9 +151,160 @@ c     instantaneous temperature
 c
 c     velocity scaling factor: Berendsen scaling
       ratio_temp = targ_temp/inst_temp
-      scalefac = dsqrt(1.d0 + (dt/tau)*(ratio_temp-1.d0))
+      scalefac = dsqrt(1.d0 + (dt/(2.d0*tau))*(ratio_temp-1.d0))
 c
-      do i = 1,nat
+      do i = 1,n
+        v(1,i) = scalefac*v(1,i)
+        v(2,i) = scalefac*v(2,i)
+        v(3,i) = scalefac*v(3,i)
+      end do
+c
+      return
+      end
+c
+c     qmd_langevin: Langevin thermostat: 
+c     Bussi and Parrinello, Phys. Rev. E (2007)
+c
+      subroutine qmd_langevin(n,dt,m,v,ekin)
+c
+      implicit none
+c
+#include "errquit.fh"
+#include "mafdecls.fh"
+#include "stdio.fh"
+#include "global.fh"
+c
+#include "qmd_common.fh"
+c
+      integer n
+      double precision dt
+      double precision m(n)    ! mass
+      double precision v(3,n)
+      double precision ekin
+c
+      integer i,j
+      double precision inst_temp
+      double precision ratio_temp
+      double precision c1,c2
+c
+      external qmd_gasdev
+      double precision qmd_gasdev
+c
+      character*32 pname
+      pname='qmd_langevin: '
+c
+c     calculate c1
+      c1 = exp(-friction*dt/2.d0)   
+      do i = 1,n
+        c2=sqrt((1.d0-c1**2)*(kb*targ_temp)/(m(i)*amu2au)) 
+        do j = 1,3
+         v(j,i) = c1*v(j,i)+c2*qmd_gasdev(idum)
+        end do
+      end do
+c
+      return
+      end
+c
+c
+c     qmd_rescale: velocity rescaling
+c
+      subroutine qmd_rescale(n,m,v,ekin,ekin2)
+c
+      implicit none
+c
+#include "errquit.fh"
+#include "mafdecls.fh"
+#include "stdio.fh"
+#include "global.fh"
+c
+#include "qmd_common.fh"
+c
+      integer n
+      double precision m(n)
+      double precision v(3,n)
+      double precision ekin
+      double precision ekin2
+c
+      integer i
+      double precision inst_temp
+      double precision ratio_temp
+      double precision scalefac
+c
+      character*32 pname
+      pname='qmd_rescale: '
+c
+c     instantaneous temperature
+      inst_temp = 2.0d0*ekin/kb/(1.d0*ndeg)
+c
+c     velocity rescaling
+c      ratio_temp = targ_temp/inst_temp
+      ratio_temp = ekin/ekin2
+      scalefac = dsqrt(ratio_temp)
+c
+      do i = 1,n
+        v(1,i) = scalefac*v(1,i)
+        v(2,i) = scalefac*v(2,i)
+        v(3,i) = scalefac*v(3,i)
+      end do
+c
+      return
+      end
+c
+c
+c     qmd_BDP: canonical velocity rescaling
+c     Bussi, Donadio, and Parrinello
+c     J. Chem. Phys. 126, 014101 (2007)
+c
+      subroutine qmd_BDP(n,dt,m,v,ekin)
+c
+      implicit none
+c
+#include "errquit.fh"
+#include "mafdecls.fh"
+#include "stdio.fh"
+#include "global.fh"
+c
+#include "qmd_common.fh"
+c
+      integer n
+      double precision dt
+      double precision m(n)
+      double precision v(3,n)
+      double precision ekin
+c
+      integer i
+      double precision scheck
+      double precision sigma
+      double precision rr, ss, cc
+      double precision scalefac
+c
+      external qmd_gasdev
+      double precision qmd_gasdev
+      external qmd_sumno
+      double precision qmd_sumno
+c
+      character*32 pname
+      pname='qmd_BDP: '
+c
+c     idum, tau, ndeg, targ_temp in common block
+c
+      if (tau.ne.0.d0) then
+        cc = exp(-dt/(2.d0*tau))
+      else
+        cc = 0.d0
+      end if
+c
+      rr = qmd_gasdev(idum)
+      sigma = dble(ndeg)*kb*targ_temp/2.d0
+      ss = qmd_sumno(ndeg-1,idum)
+      scheck = rr + dsqrt(cc*dble(ndeg)*ekin/((1.d0-cc)*sigma))
+c
+      scalefac = cc + (sigma/(dble(ndeg)*ekin))*(1.d0-cc)*(ss+rr**2) 
+     $          + 2.d0*rr*dsqrt(cc*(1.d0-cc)*sigma/(dble(ndeg)*ekin))
+c
+      scalefac = dsign(dsqrt(scalefac),scheck)
+c
+      do i = 1,n
         v(1,i) = scalefac*v(1,i)
         v(2,i) = scalefac*v(2,i)
         v(3,i) = scalefac*v(3,i)
diff --git a/src/qmd/qmd_write_rst.F b/src/qmd/qmd_write_rst.F
index 13b4d58cdb..8a1d5c28bb 100644
--- a/src/qmd/qmd_write_rst.F
+++ b/src/qmd/qmd_write_rst.F
@@ -2,7 +2,7 @@ c
 c     write restart
 c
       subroutine qmd_write_rst(rstfilename,rstfileunit,
-     &  nat,iframe,esys,ekin,geom)
+     &  n,iframe,ekin,geom)
 c
       implicit none
 c
@@ -17,14 +17,14 @@ c
 c
       character*(*) rstfilename
       integer rstfileunit
-      integer nat
+      integer n
       integer iframe
-      double precision esys,ekin,etotal
+      double precision ekin
       integer geom
 c
       integer i
-      double precision r(3,nat)
-      double precision v(3,nat)
+      double precision r(3,n)
+      double precision v(3,n)
 c
       character*32 pname
 c
@@ -35,29 +35,27 @@ c     open file
      &        position='rewind',status='unknown')
 c
 c     write number of atoms, frame number and energy
-      write(rstfileunit,110) nat
-      etotal = esys + ekin
-      write(rstfileunit,111) iframe,etotal,esys,ekin
-      write(rstfileunit,112) targ_temp,final_temp
+      write(rstfileunit,110) n
+      write(rstfileunit,111) iframe,ekin
       if (.not.geom_cart_coords_get(geom,r))
      $  call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
       if (.not.geom_vel_get(geom,v))
      $   call errquit(pname//'geom_vel_set',0,GEOM_ERR)
 c
       write(rstfileunit,*) "Coordinates"
-      do i=1,nat
+      do i=1,n
          write(rstfileunit,113) r(1,i),r(2,i),r(3,i)
       end do
 c
       write(rstfileunit,*) "Velocities"
-      do i=1,nat
+      do i=1,n
          write(rstfileunit,113) v(1,i),v(2,i),v(3,i)
       end do
 c
       if (rstfileunit.ne.0) close(unit=rstfileunit)
 c
   110 format(I5)
-  111 format(I5,3f22.12)
+  111 format(I10,3f22.12)
   112 format(2f22.12)
   113 format(3f22.12)
 c
diff --git a/src/qmd/qmd_write_trj.F b/src/qmd/qmd_write_trj.F
index 70f022092d..3a9987e294 100644
--- a/src/qmd/qmd_write_trj.F
+++ b/src/qmd/qmd_write_trj.F
@@ -2,7 +2,7 @@ c
 c     qmd_write_trj: write trajectory: xyz file in angstroms
 c
       subroutine qmd_write_trj(xyzfilename,xyzfileunit,
-     &           mode_append,nat,iframe,etotal,geom)
+     &           mode_append,n,iframe,etotal,geom)
 c
       implicit none
 c
@@ -18,10 +18,15 @@ c
       character*(*) xyzfilename
       integer xyzfileunit
       logical mode_append
-      integer nat
+      integer i, j
+      integer n
       integer iframe
       double precision etotal
       integer geom
+      character(16) t(n)
+      double precision r(3,n)
+      double precision v(3,n)
+      double precision q(n)
 c
       character*32 pname
 c
@@ -37,16 +42,25 @@ c     open file
       end if
 c
 c     write number of atoms, frame number and energy
-      write(xyzfileunit,110) nat
+      write(xyzfileunit,110) n
       write(xyzfileunit,111) iframe,etotal
 c
-      if (.not.geom_print_xyzq(geom, xyzfileunit)) 
-     &  call errquit(pname//'geom write failed',0,GEOM_ERR)
+c     get number of centers, tags, coordinates, and charges
+      if (.not.geom_cart_get(geom,i,t,r,q))
+     $  call errquit(pname//'geom_cart_coords_get',0,GEOM_ERR)
+c     get velocities
+      if (.not.geom_vel_get(geom,v))
+     $   call errquit(pname//'geom_vel_set',0,GEOM_ERR)
+c
+      do i=1,n
+       write(xyzfileunit,4) t(i),r(:,i)*au2ang,v(:,i)
+      end do
 c
       if (xyzfileunit.ne.0) close(unit=xyzfileunit)
 c
-  110 format(I5)
-  111 format(I5,f22.12)
+ 4    format(1x,a12,3f15.8,2x,3f15.8)
+ 110  format(I5)
+ 111  format(I10,f22.12)
 c
       return
       end