Accessing electron density from dplot via python

property
  esp
  grid output esp.cube
end
.....
python
	#Carry out some rotations
        ....
	# Update file name
	rtdb_put('prop:grid:output', "0_000" +'.cube')
	print task_property('dft')
end

dplot
...
end

....  
python
	#Carry out some rotations
	....
        # Update file name
        rtdb_put('dplot:File_Out', "0_000" +'.density.cube')
	print dplot
end

--- /tmp/Nwchem-6.5.revision26243-src.2014-09-10/src/property/hnd_elpmap.F      2014-09-10 19:08:56.000000000 +0100
+++ ./property/hnd_elpmap.F     2015-04-27 13:39:43.997092910 +0100
@@ -24,6 +24,7 @@ c
       character*16 element, at_tag
       integer iat, atn, nat, i
       integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt, l_epot, k_epot
+      integer l_eden, k_eden
       integer nefc, l_efcc, k_efcc, l_efcz, k_efcz
       integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2)
       character*3 scftyp
@@ -160,12 +161,21 @@ c     total electric field array (M.V.)
 c     --------------------------------
       if (.not. ma_push_get(mt_dbl,nprp,'tot epot',l_tepot,k_tepot))
      &    call errquit('hnd_elfmap: ma failed',911,MA_ERR)
+
+      if (.not. ma_push_get(mt_dbl,nprp,'eden',l_eden,k_eden))
+     &    call errquit('hnd_elpmap: eden ma failed',911,MA_ERR)
+
+
 cc
 c     ----- calculate electronic contribution at all points -----
 c
       call hnd_elfcon(basis,geom,g_dens(ndens),dbl_mb(i_prp_c),nprp,
      &                dbl_mb(k_epot),0)

+c     ----- calculate electron density                     ------
+      call hnd_elfcon(basis,geom,g_dens(ndens),dbl_mb(i_prp_c),nprp,
+     &                dbl_mb(k_eden),-1)
+
       dobq = .false.
       if(geom_extbq_on()) then
         dobq = .true.
@@ -311,6 +321,14 @@ c
         call util_print_centered(luout,
      >    "writing esp potential to "//cube_name,.true.,.true.)
         call prop_grid_write_cube(geom,nprp,dbl_mb(k_tepot),cube_name)
+
+        if (.not. rtdb_cget(rtdb, "prop:grid:output",1,cube_name))
+     >     call util_file_prefix("eld.cube",cube_name)
+
+        cube_name = trim(cube_name)//'_eld.cube'
+        call util_print_centered(luout,
+     >    "writing electron density to "//cube_name,.true.,.true.)
+        call prop_grid_write_cube(geom,nprp,dbl_mb(k_eden),cube_name)
       end if

 c
@@ -325,10 +343,13 @@ c     ----------------------------------
 cc
 c     ------- Deallocate MA memory ------
 c
+      if (.not.ma_pop_stack(l_eden)) call errquit
+     &   ('hnd_elpmap, ma_pop_stack of l_eden failed',911,MA_ERR)
       if (.not.ma_pop_stack(l_tepot)) call errquit
      &   ('hnd_elpmap, ma_pop_stack of l_tepot failed',911,MA_ERR)
       if (.not.ma_pop_stack(l_epot)) call errquit

dplot
  GAUSSIAN
  TITLE Density
   LimitXYZ
   -4.611 4.718 107
   -5.058 4.295 108
   -4.042 4.042 93
  spin total
  output density.cube
end