OK, it would appear that I may have some form of a solution in a patch to hnd_elpmap. Does the following approach seem reasonable?
--- /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
|