diff --git a/src/geom/geom.F b/src/geom/geom.F
index 282989e107..8062854202 100644
--- a/src/geom/geom.F
+++ b/src/geom/geom.F
@@ -1814,6 +1814,7 @@ c
       end do
 c
       end
+
 c
 C> \brief Extract only the coordinates and charges from a geometry
 C> instance
@@ -2762,8 +2763,11 @@ c
          call util_print_centered(LuOut,'Velocities',40,.true.)
          write(LuOut,3) 
          do icent = 1, ncenter(geom)
+            do i = 1, 3
+               tmp(i) = velocities(i,icent,geom)/scale ! Scale units as necessary
+            end do
             write(LuOut,4) icent, tags(icent,geom), charge(icent,geom),
-     $           (velocities(i,icent,geom),i=1,3)
+     $           (tmp(i),i=1,3)
          end do
       endif
 c
diff --git a/src/geom/geom_input.F b/src/geom/geom_input.F
index 364f28208d..edff6b9a0c 100644
--- a/src/geom/geom_input.F
+++ b/src/geom/geom_input.F
@@ -677,6 +677,9 @@ c
       if (.not. geom_cart_get(geom, ncenter, tags, coords, charge))
      $     call errquit('geom_input: geom_cart_get failed', 0,
      &       GEOM_ERR)
+      if (.not. geom_vel_get(geom, velocities))
+     $     call errquit('geom_input: geom_vel_get failed', 0,
+     &       GEOM_ERR)
 
 c
 c     If doing user zmatrix then remember values of constant variables
diff --git a/src/nwpw/nwpwlib/ion/ion.F b/src/nwpw/nwpwlib/ion/ion.F
index c2bd622702..f448546cae 100644
--- a/src/nwpw/nwpwlib/ion/ion.F
+++ b/src/nwpw/nwpwlib/ion/ion.F
@@ -3,9 +3,9 @@
 *
 
 *     ***************************
-*     *				*
-*     *		ion_end 	*
-*     *				*
+*     *                         *
+*     *        ion_end          *
+*     *                         *
 *     ***************************
       subroutine ion_end()
       implicit none
@@ -54,9 +54,9 @@
       end
 
 *     ***************************
-*     *				*
-*     *		ion_write	*
-*     *				*
+*     *                         *
+*     *        ion_write        *
+*     *                         *
 *     ***************************
       subroutine ion_write(rtdb)
       implicit none
@@ -164,9 +164,9 @@
       end
 
 *     ***************************
-*     *				*
-*     *		ion_destroy 	*
-*     *				*
+*     *                         *
+*     *        ion_destroy      *
+*     *                         *
 *     ***************************
       subroutine ion_destroy()
       implicit none
@@ -187,9 +187,9 @@
     
 
 *     ***************************
-*     *				*
-*     *		ion_init	*
-*     *				*
+*     *                         *
+*     *        ion_init         *
+*     *                         *
 *     ***************************
 *
 *        This routine reads in the ion data structure from the
@@ -1991,9 +1991,9 @@ c      end if
 
 
 *     ***************************
-*     *				*
-*     *		ion_vion	*
-*     *				*
+*     *                         *
+*     *          ion_vion       *
+*     *                         *
 *     ***************************
 
       double precision function ion_vion(i,j)
diff --git a/src/symmetry/sym_nwc.F b/src/symmetry/sym_nwc.F
index ee1c8e9453..f68cc74590 100644
--- a/src/symmetry/sym_nwc.F
+++ b/src/symmetry/sym_nwc.F
@@ -78,10 +78,17 @@ c
      &     latags, iatags)
       LResult = MA_Push_Get(MT_Dbl,nata*(nops+1)*3,'scratch atoms',
      &     lcoord, icoord) .and. lresult
+      LResult = MA_Push_Get(MT_Dbl,nata*(nops+1)*3,'scratch velocities',
+     &     lvel, ivel) .and. lresult
       Lresult = MA_Push_Get(MT_Int,nata,'distinct atoms',ldatom,idatom)
      $     .and. lresult
       if (.not. lresult) call errquit('sym_nwc: dosymops ma failed', 0,
      &       MA_ERR)
+      do ii=1,nata
+          dbl_mb(ivel+3*(ii-1))   = velocities(1,ii,geom)
+          dbl_mb(ivel+3*(ii-1)+1) = velocities(2,ii,geom)
+          dbl_mb(ivel+3*(ii-1)+2) = velocities(3,ii,geom)
+      end do
 c
       call dosymops(sym_ops(1,1,geom),nops,coords(1,1,geom),
      $     tags(1,geom),nata,
@@ -89,6 +96,7 @@ c
      $     cdist,cang,Int_MB(idatom), 
      $     charge(1,geom), charge_new, geom_mass(1,geom),
      $     mass_new, geom_invnucexp(1,geom),invnucexp_new,threquiv)
+
 c
       if(nareal.gt.max_cent) then
         write(6,*)' Too many atoms, increase nw_max_atom in ',
@@ -139,7 +147,7 @@ c     remain in icoord.
 c     Note that the memory for the map (lmscr) is now freed in geom_destroy
 c
       call sym_put_in_geom(geom, nata, nareal,
-     $     Byte_MB(iatags), int_mb(idatom), dbl_mb(icoord),
+     $     Byte_MB(iatags), int_mb(idatom),dbl_mb(icoord),dbl_mb(ivel),
      $     lmscr, imscr, nops, charge_new, mass_new, invnucexp_new)
 c
 c ... sloppy way for the moment to fix units for 2-d case
@@ -221,6 +229,7 @@ c
 c--> final clean up of memory: gensym and dosymops activities
 c
       LResult = MA_Pop_Stack(ldatom)
+      LResult = MA_Pop_Stack(lvel)
       LResult = MA_Pop_Stack(lcoord)
       LResult = MA_Pop_Stack(latags)
 c
diff --git a/src/symmetry/sym_put_geom.F b/src/symmetry/sym_put_geom.F
index 30e7dcb9a1..642f41cdd0 100644
--- a/src/symmetry/sym_put_geom.F
+++ b/src/symmetry/sym_put_geom.F
@@ -1,5 +1,6 @@
       subroutine sym_put_in_geom(geom, nat_uniq, nat_new,
-     $     tags_new, uniq_cent, coords_new, l_map, k_map, nops,
+     $     tags_new, uniq_cent, coords_new, velocities_new,
+     $     l_map, k_map, nops,
      $     charge_new, mass_new, invnucexp_new)
 C$Id$
       implicit none
@@ -14,6 +15,7 @@ C$Id$
       character*16 tags_new(nat_new)
       integer uniq_cent(nat_uniq)
       double precision coords_new(3,nat_new)
+      double precision velocities_new(3,nat_new)
       double precision charge_new(max_cent)
       double precision mass_new(max_cent)
       double precision invnucexp_new(max_cent)
@@ -34,6 +36,9 @@ c
       enddo
       call dgemm('n', 'n', 3, nat_new, 3, 1.0d0, amatrix(1,1,geom), 3,
      $     coords_new, 3, 0.0d0, coords(1,1,geom), 3)
+      call dgemm('n', 'n', 3, nat_new, 3, 1.0d0, amatrix(1,1,geom), 3,
+     $     velocities_new, 3, 0.0d0, velocities(1,1,geom), 3)
+c
 c
       ncenter_unique(geom) = nat_uniq
       ncenter(geom) = nat_new