diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F
index ec6fa43c42..bf6b74d0fb 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F
@@ -616,7 +616,7 @@ c
         call tddft_grad_compute_dao(ipol,nao,nocc,g_mo,g_d)
 c DEBUG
 c        call ga_print(g_d)
-        if (tddft_grad_util_print('tddft grad g',print_debug)) then
+        if (util_print('tddft grad g',print_debug)) then
           if (oroot) write(LuOut,*)'DEBUG: '//pname//'D'
           call tddft_grad_print_array(ipol,1,g_d,dble(ipol))
         endif
@@ -640,7 +640,7 @@ c
         enddo
 c DEBUG
 c        call ga_print(g_p_ao)
-        if (tddft_grad_util_print('tddft grad g',print_debug)) then
+        if (util_print('tddft grad g',print_debug)) then
           if (oroot) write(LuOut,*)'DEBUG: '//pname//'P'
 c          call tddft_grad_print_array(ipol,1,g_p_ao,dble(ipol))
           call tddft_grad_print_array(ipol,1,g_d(3),dble(ipol))
@@ -798,7 +798,7 @@ c
 c Daniel (11-30-12): Print X+Y and X-Y in the AO basis
 c DEBUG
 c        call ga_print(g_d(5))
-        if (tddft_grad_util_print('tddft grad g',print_debug)) then
+        if (util_print('tddft grad g',print_debug)) then
           if (oroot) write(LuOut,*)'DEBUG: '//pname//'(X+Y)'
           call tddft_grad_print_array(ipol,1,g_d(5),dble(ipol))
         endif
@@ -818,7 +818,7 @@ c Daniel (11-30-12): Print W (note, this contains the ground state
 c density matrix contribution as well).
 c DEBUG
 c        call ga_print(g_wp_ao)
-        if (tddft_grad_util_print('tddft grad g',print_debug)) then
+        if (util_print('tddft grad g',print_debug)) then
           if (oroot) write(LuOut,*)'DEBUG: '//pname//'W'
           call tddft_grad_print_array(1,1,g_wp_ao,dble(ipol))
         endif
@@ -1140,7 +1140,7 @@ c Daniel (4-8-13): Need conditional statement here for CAM functionals
 c
 c       Print results if so required
 c
-        if (tddft_grad_util_print('tddft grad terms',print_high)) then
+        if (util_print('tddft grad terms',print_high)) then
           if (oroot) then
             write(LuOut,5)ir
     5       format(' The Excited State Energy Gradient by Terms',
diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_p.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_p.F
index 15aef30c36..efe230cc89 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_p.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_p.F
@@ -102,7 +102,7 @@ c
         call nga_scale_patch(g_tp(ip),blo,bhi,0.5d0)
       enddo
 c
-      if (tddft_grad_util_print(pname,print_debug)) then
+      if (util_print(pname,print_debug)) then
         oroot = ga_nodeid().eq.0
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'P'
         call tddft_grad_print_array(ipol,nroots,g_tp,dble(ipol))
diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_r.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_r.F
index 14daed1c75..a2a80a1a41 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_r.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_r.F
@@ -120,17 +120,21 @@ c
 c
 c Daniel (2-8-13): This is a local TDA variable, so that we don't need
 c to change the global tda variable.
-      logical tdaloc
+      logical tdaloc, tdaloc_rtdb
+      logical oprint
+      integer g_r_x(2),g_r_org(2)
 c
       double precision Exc(2)    ! Exchange-correlation energy
 c
       pname="tddft_grad_compute_r: "
       iwhich = 0  ! call to tddft_nga_cont()
+      write(6,*) ' 0 tda ',tda
 c
 c     1. Do the T contributions
 c
 c     1.1. Create global array to hold Tuv
 c
+      oprint=util_print('tddft_grad_r',print_debug)
       idim(1) = nroot*ipol
       idim(2) = nao
       idim(3) = nao
@@ -180,6 +184,7 @@ c fock_xc:triplet here for restricted triplet calculations to work.
 c
       call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_tuv,g_apbt,g_ambt,
      +nao,ipol,tol2e,tda,oskel,kfac,lhashf,.false.,nroot,iwhich,lhascd)
+      write(6,*) ' 2 tda ',tda
 c Daniel (2-26-13): Reset fock_xc:triplet here for restricted triplet 
 c calculations to work.
       if (otriplet) then
@@ -217,11 +222,11 @@ c
 c
 c     1. Done T contributions
 cDEBUG
-c     if (tddft_grad_util_print('tddft grad r',print_debug)) then
-c       oroot = ga_nodeid().eq.0
-c       if (oroot) write(LuOut,*)'DEBUG A: tddft_grad_comp_r: R'
-c       call tddft_grad_print_array(ipol,nroot,g_r,dble(ipol))
-c     endif
+      if (oprint) then
+       oroot = ga_nodeid().eq.0
+       if (oroot) write(LuOut,*)'DEBUG A: tddft_grad_comp_r: R_T'
+       call tddft_grad_print_array(ipol,nroot,g_r,dble(ipol))
+      endif
 cDEBUG
 c
 c     2. Do the (X+Y) and (X-Y) contributions
@@ -304,6 +309,8 @@ c the global one.
           tdaloc = .false. ! For TDDFT calculations
         endif
       endif
+        if (rtdb_get(rtdb,'tddft:tdaloc',mt_log,1,tdaloc_rtdb))
+     +    tdaloc=tdaloc_rtdb
       call tddft_nga_cont(rtdb,ihdl_geom,ihdl_bfao,g_x,g_apbx,g_ambx,
      +nao,ipol,tol2e,tdaloc,oskel,kfac,lhashf,otriplet,nroot,iwhich,
      +lhascd)
@@ -349,6 +356,21 @@ c     2.10 Transform (A+B)(X+Y) to MO basis occupied-occupied block only
 c
       call tddft_grad_trans_ao2mo(ipol,nao,nfc,naoc,nocc,nav,nfv,
      +     nroot,1.0d0,0.0d0,"ij",g_mo,g_apbx,g_hij,"ij")
+      if(oprint)then
+c     extra ga to keep track of X
+         if (.not. ga_duplicate(g_r(1), g_r_x(1), 'gtmp2'))
+     c        call errquit('tdgcr: gadupl failed',0,0)
+         call ga_zero(g_r_x(1))
+         g_r_org(1)=g_r(1)
+         g_r(1)=g_r_x(1)
+         if(ipol.eq.2) then
+         if (.not. ga_duplicate(g_r(2), g_r_x(2), 'gtmp2'))
+     c           call errquit('tdgcr: gadupl failed',0,0)
+         call ga_zero(g_r_x(2))
+         g_r_org(2)=g_r(2)
+         g_r(2)=g_r_x(2)
+         endif
+      endif
 c
 c     2.11 Add -sum_j (X+Y)ja [(A+B)(X+Y)ji] to Ria
 c
@@ -546,11 +568,21 @@ c
 c     2. Done (X+Y) and (X-Y) contributions
 c
 cDEBUG
-c     if (tddft_grad_util_print('tddft grad r',print_debug)) then
-c       oroot = ga_nodeid().eq.0
-c       if (oroot) write(LuOut,*)'DEBUG B: tddft_grad_comp_r: R'
-c       call tddft_grad_print_array(ipol,nroot,g_r,dble(ipol))
-c     endif
+      if (oprint) then
+       oroot = ga_nodeid().eq.0
+       if (oroot) write(LuOut,*)'DEBUG B: tddft_grad_comp_r: R_x'
+       call tddft_grad_print_array(ipol,nroot,g_r_x,dble(ipol))
+       call ga_add(1d0,g_r_x(1),1d0,g_r_org(1),g_r_org(1))
+       g_r(1)=g_r_org(1)
+       if (.not.ga_destroy(g_r_x(1))) 
+     +     call errquit(pname//'failed to destroy grx',1,GA_ERR)
+       if(ipol.eq.2) then
+       call ga_add(1d0,g_r_x(2),1d0,g_r_org(2),g_r_org(2))
+       g_r(2)=g_r_org(2)
+       if (.not.ga_destroy(g_r_x(2))) 
+     +     call errquit(pname//'failed to destroy grx',2,GA_ERR)
+       endif
+      endif
 cDEBUG
 c
 c     3. Do the Gxc contributions
@@ -664,14 +696,14 @@ c
         enddo
         call ga_sync()
 cDEBUG
-c     if (tddft_grad_util_print('tddft grad r',print_debug)) then
-c       oroot = ga_nodeid().eq.0
-c       if (oroot) write(LuOut,*)'DEBUG C: tddft_grad_comp_r: den'
-c       call tddft_grad_print_array(ipol,nroot,int_mb(k_den2),
-c    +                              dble(ipol))
-c       call tddft_grad_print_array(ipol,nroot,int_mb(k_den2+1),
-c    +                              dble(ipol))
-c     endif
+        if (oprint) then
+           oroot = ga_nodeid().eq.0
+           if (oroot) write(LuOut,*)'DEBUG C: tddft_grad_comp_r: den'
+       call tddft_grad_print_array(ipol,nroot,int_mb(k_den2),
+     +          dble(ipol))
+       call tddft_grad_print_array(ipol,nroot,int_mb(k_den2+1),
+     +      dble(ipol))
+      endif
 cDEBUG
 c Daniel (1-10-13): For a single restricted root ngxc = 1 and there are 2
 c density matrices fed to fock_xc, the GS density matrix D and the
@@ -689,7 +721,7 @@ c fock_xc.
      +               int_mb(k_gxc), Exc, ipol, .false.)
         call ga_sync()
 c DEBUG
-      if (tddft_grad_util_print('tddft grad r',print_debug)) then
+      if (oprint) then
         oroot = ga_nodeid().eq.0
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'gxc'
         call tddft_grad_print_array(ipol,nroot,int_mb(k_gxc),
@@ -734,9 +766,9 @@ c
       endif  ! xc_gotxc()
 c
 c     Debug print
-      if (tddft_grad_util_print('tddft grad r',print_debug)) then
+      if (oprint) then
         oroot = ga_nodeid().eq.0
-        if (oroot) write(LuOut,*)'DEBUG: '//pname//'R'
+        if (oroot) write(LuOut,*)'DEBUG: '//pname//'R_all'
         call tddft_grad_print_array(ipol,nroot,g_r,dble(ipol))
       endif
 c
diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_t.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_t.F
index 756a133b92..e89ebd5dee 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_t.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_t.F
@@ -147,7 +147,7 @@ c (X-Y)ias(X-Y)ibs
 c
       enddo ! ip = 1, ipol
 c
-      if (tddft_grad_util_print('tddft grad t',print_debug)) then
+      if (util_print('tddft grad t',print_debug)) then
         oroot = ga_nodeid().eq.0
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'T'
         call tddft_grad_print_array(ipol,nroot,g_tp,dble(ipol))
diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_w.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_w.F
index c411e13e82..b43f014203 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_w.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_w.F
@@ -599,7 +599,7 @@ c Daniel (2-12-13): We need to set triplet here for fock_xc.
      +               int_mb(k_gxc), Exc, ipol, .false.)
         call ga_sync()
 c DEBUG
-      if (tddft_grad_util_print('tddft grad w',print_debug)) then
+      if (util_print('tddft grad w',print_debug)) then
         oroot = ga_nodeid().eq.0
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'gxc'
         call tddft_grad_print_array(ipol,nroots,int_mb(k_gxc),
@@ -1190,7 +1190,7 @@ c
       enddo
       call ga_sync()
 c
-      if (tddft_grad_util_print('tddft grad w',print_debug)) then
+      if (util_print('tddft grad w',print_debug)) then
         oroot = ga_nodeid().eq.0
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'W'
         call tddft_grad_print_array(ipol,nroots,g_w,dble(ipol))
diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_solve.F b/src/nwdft/lr_tddft_grad/tddft_grad_solve.F
index 258d8c8b78..9313f66715 100644
--- a/src/nwdft/lr_tddft_grad/tddft_grad_solve.F
+++ b/src/nwdft/lr_tddft_grad/tddft_grad_solve.F
@@ -254,7 +254,7 @@ c
       if (.not.ga_destroy(g_x)) 
      +   call errquit(pname//'failed to destroy g_x',0, GA_ERR)
 c
-      if (tddft_grad_util_print('tddft grad z',print_debug)) then
+      if (util_print('tddft grad z',print_debug)) then
         if (oroot) write(LuOut,*)'DEBUG: '//pname//'Z'
         call tddft_grad_print_array(ipol,nroots,g_z,dble(ipol))
       endif
diff --git a/src/nwdft/lr_tddft_grad/tddft_gradient.F b/src/nwdft/lr_tddft_grad/tddft_gradient.F
index 928688ab70..2794d8211b 100644
--- a/src/nwdft/lr_tddft_grad/tddft_gradient.F
+++ b/src/nwdft/lr_tddft_grad/tddft_gradient.F
@@ -168,6 +168,9 @@ c     Get the start times
 c
 c     Init pstat
       call dft_pstat_init(rtdb)
+c      
+      call util_print_push
+      call util_print_rtdb_load(rtdb,'tddft')
 c
 c     Print header
       if (oroot) call tddft_grad_print_header(luout)
@@ -225,7 +228,7 @@ c     Print calculation information
 c
 c     Print start times
       if (oroot) then
-        if (tddft_grad_util_print('tddft grad time',print_default)) then
+        if (util_print('tddft grad time',print_default)) then
           write(luout,*)
           write(luout,"(' Start at time  cpu:',f11.1,'s  wall:',f11.1,
      +                  's')")-tddft_grad_cpu,-tddft_grad_wall
@@ -435,7 +438,7 @@ c     Print end times and wrap up
       tddft_grad_cpu  = tddft_grad_cpu  + util_cpusec()
       tddft_grad_wall = tddft_grad_wall + util_wallsec()
       if (oroot) then
-        if (tddft_grad_util_print('tddft grad time',print_default)) then
+        if (util_print('tddft grad time',print_default)) then
           write(LuOut,*)
           write(LuOut,"(' TDDFT Gradient time  cpu:',f11.1,'s  wall:',
      +                  f11.1,'s')")tddft_grad_cpu,tddft_grad_wall
@@ -496,6 +499,7 @@ c
 c <-- MN solvation models 
 c
       call dft_pstat_print
+      call util_print_pop
       return
       end
 c $Id$