new geom keyword store_symrot to print gradients with the original geometry before autosym

src/geom/geom_input.F

@@ -32,6 +32,7 @@ c
       logical ozmatrix
       logical oautoz
       logical oautosym
+      logical store_symrot
       logical ocenter
       logical found_cart
       logical ofinite, oatomfinite
@@ -92,8 +93,10 @@ c
       logical found_ang2au_conv
       logical found_au2ang_conv
       double precision ang2au, au2ang
+      double precision rot(3,3),rot_inv(3,3)
       integer istart,iend,irate
       integer atn, isys, fcent
+      integer l_xyzold,k_xyzold
       character*16 format, nucmodel
       logical geom_zmtmak
       external geom_zmtmak
@@ -122,7 +125,7 @@ c
 c     
 c     geometry [<name>] [units <units>] [print] [noprint] [bqbq] \
 c     [angstrom_to_au <real>] [au_to_angstrom <real>] \
-c     [(no)autoz] [nucleus <model>]
+c     [(no)autoz] [nucleus <model>] [store_symrot]
 c     
 c     subdirectives
 c     
@@ -161,6 +164,7 @@ c
       ecce_periodic     = .false.
       sysset            = .false.
       symset            = .false.
+      store_symrot      = .false.
 c
       call dfill(3*max_center, 0.0d0, velocities, 1)
 c
@@ -257,6 +261,8 @@ c ------ Store ofinite to be used by HFine finite calc ---FA-03-21-11-- START
      $      call errquit('geo ofinite: rtdb_put failed',555, RTDB_ERR)
 c ------ Store ofinite to be used by HFine finite calc ---FA-03-21-11-- END
             end if
+         else if (inp_compare(.false.,'store_symrot',field)) then
+            store_symrot = .true.
          else
             if (name .ne. ' ') call errquit
      $           ('geom_input: geometry [<name>] [units <units>]', 0,
@@ -654,10 +660,25 @@ c
             write(luout,*)
          endif
 c     
+         if(store_symrot) then
+            if(.not. ma_push_get(mt_dbl,3*ncenter,'xyzold',
+     L        l_xyzold,k_xyzold)) call
+     $           errquit('geom_input: ma error', ncenter, MA_ERR)
+            call
+     D           dcopy(ncenter*3, coords, 1, dbl_mb(k_xyzold), 1)
+         endif
          call geom_getsym0(rtdb,geom,oautosym,
      ,     ncenter,isys,
      ,     name,tags,
      ,     threquiv,scale,coords,charge,velocities,atomct)
+         if(store_symrot) then
+            call geom_findrot(ncenter,dbl_mb(k_xyzold),coords,rot)
+            call geom_m3inv(rot,rot_inv)
+            if (.not. rtdb_put(rtdb, 'geom:symrot', MT_DBL, 9, rot_inv)) 
+     E           call errquit("geom_input: symrot",0,0)
+            if (.not. ma_pop_stack(l_xyzold))
+     &           call errquit('geom_input:cant pop stack',0,0)
+         endif
       endif                     ! End of else clause for if(oadjust)
 
 c

src/gradients/grad_force.F

@@ -168,14 +168,17 @@ c
       integer ifocc
       logical oskel, omp2, odft, ocdfit, status,frac_occ
       double precision xfac(numfunc), jfac, kfac
+      double precision rot(3,3),xrot,yrot,zrot
 c
       logical has_frac_occ
       external has_frac_occ
 c
+      double precision fol(3),f1,f2,f3
       logical lcdft  ! Qin
       integer ipop
       logical out1
       logical zerof
+      logical have_symrot
 c
 cc AJL/Begin/SPIN ECPs
       integer ecp_channels
@@ -233,7 +236,9 @@ c
           do_cosmo_smd = .false.
         endif
       endif
-
+c     sym rot matrix
+      have_symrot=rtdb_get(rtdb, 'geom:symrot', MT_DBL, 9, rot)
+     
 C     get SCF MO vectors for density
 
       odft = .false.
@@ -1234,8 +1239,26 @@ C     default: print the total forces
          do 30, i=1, nat
             if (.not. geom_cent_get(geom, i, tag, crd, q)) call errquit
      $           ('gradients: geometry corrupt?',0, GEOM_ERR)
+            if(have_symrot) then
+         xrot = crd(1)*rot(1,1) + crd(2)*rot(2,1) + crd(3)*rot(3,1)
+         yrot = crd(1)*rot(1,2) + crd(2)*rot(2,2) + crd(3)*rot(3,2)
+         zrot = crd(1)*rot(1,3) + crd(2)*rot(2,3) + crd(3)*rot(3,3)
+         crd(1)=xrot
+         crd(2)=yrot
+         crd(3)=zrot
+         fol(1)=dbl_mb(k_force+3*(i-1))
+         fol(2)=dbl_mb(k_force+3*(i-1)+1)
+         fol(3)=dbl_mb(k_force+3*(i-1)+2)
+         f1 = fol(1)*rot(1,1) + fol(2)*rot(2,1) + fol(3)*rot(3,1)
+         f2 = fol(1)*rot(1,2) + fol(2)*rot(2,2) + fol(3)*rot(3,2)
+         f3 = fol(1)*rot(1,3) + fol(2)*rot(2,3) + fol(3)*rot(3,3)
+         else
+            f1=dbl_mb(k_force+3*(i-1))
+            f2=dbl_mb(k_force+3*(i-1)+1)
+            f3=dbl_mb(k_force+3*(i-1)+2)
+         endif
             write(luout,2000) i, tag,(crd(j),j=1,3),
-     $           (dbl_mb(k_force+3*(i-1)+j),j=0,2)
+     $           f1,f2,f3
  30      continue
          write(luout,*)
  1000    format(/,/,25X,A,' ENERGY GRADIENTS',/,/,4X,'atom',15X,