Accessing electron density from dplot via python


Clicked A Few Times
Dear All,

Within an embedded python script in a NWChem 6.5 calculation, I am rotating a molecule and obtaining
its ESP at these different orientations:

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

However, I would also like to obtain an electron density at these conformation as well,
and following the pattern above, it seems that dplot should be able to do this:

dplot
...
end

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

However, this does not work; I understand dplot is special in some way. Is there anyway I can access the electron density within the python script in the same way that I have obtained the ESP above?

Clicked A Few Times
Looking more into the source, would it be possible to use the ELECTRONDENSITY property, i.e. the call to hnd_eldmap, and add cube writing functionality to this routine in a similar way hnd_elfmap and hnd_elpmap this ability?

Clicked A Few Times
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


Forum Vet
It seems reasonable to me.
Does it seem to produce meaningful results?

Clicked A Few Times
Thank you for inspecting the patch. Yes, comparing this density to the output from another QM code seems OK, as well as to the output from a dplot calculation (which is not supported via the python interface):

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

Visual inspection via VMD also seems OK.

A little background on why I'm doing this. I'm essentially carrying out multiple MEP calculations on the same molecule at different orientations to average discretization effects from the projection onto a grid. I am using the MEP at a specific isosurface at which the electron density has a value of 0.002 (hence the need for these two cubes). This MEP isosurface is then used for a "coarse grained" abstraction reducing the molecule down to a set of Surface Site Interaction Points (SSIPs). More about this approach can be found here: http://dx.doi.org/10.1039/C3SC22124E

I am thinking perhaps of perhaps submitting a patch so that NWChem code so that for the property esp, one can pass an option to only calculate the MEP values at a specific electron density isovalue (within a tolerance), but this may require a bit more work than the patch above to write just the electron density.

As an aside, the python interface is extremely useful for bespoke and non-standard calculations; I really appreciate this functionality.

Forum Vet
I have just checked this patch into the NWChem repository.
Could I get your name and affiliation so that the contribution does not remain anonymous?

Clicked A Few Times
Dear Edoapra,

That's great; thank you for that. I'm Mark Williamson, Cambridge University, UK mw529@cam.ac.uk

http://www-hunter.ch.cam.ac.uk/index.php?page=people

Forum Vet
Thanks for your contribution

Clicked A Few Times
Dear Edoapra,

As a follow up to this, I have implemented a few more grid related fixes and an extra option to specify the step size of the grid spacing as an input, as opposed to its internal hard-coded value of 0.2 A. When this is invoked, the cap of 50 points per grid direction is removed.

In addition, I have added a functionality so that only a subset of electrostatic potential grid points that lie on a user specified electron density surface (within a specified tolerance) are calculated. It would be really great if these could make it back to the original code.

There are quite a few changes, sprinkled over multiple files, and an extra file. I'm wondering if there is a better way to apply these rather than pasting multiple patches to this forum. In addition, I have a few queries about my changes, hence the patches may go through a few iterations. I notice Jeff has a git repository (https://github.com/jeffhammond/nwchem) that is tracking original SVN repo at https://svn.pnl.gov/svn/nwchem . In theory, I could fork Jeff's github repo at its current HEAD and then issue a pull request, however, I am not sure how easy it would then be to back propagate this to the SVN repository.

What action do you recommend?

Thanks,

Mark

Gets Around
Edo:

I'm going to help Mark with this, but I don't plan to use git-svn to commit changes to the PNNL repo.

Clicked A Few Times
Is there any update on this?

I have a pull request (but with a few questions) at:
https://github.com/jeffhammond/nwchem/pull/1

If there is anything more I can do, please let me know.

Thanks again,

Mark


Forum >> NWChem's corner >> General Topics