Accessing electron density from dplot via python

From NWChem

Viewed 303 times, With a total of 7 Posts
Jump to: navigation, search
Forum >> NWChem's corner >> General Topics

Clicked A Few Times
Threads 1
Posts 5
10:38:04 AM PDT - Mon, Apr 20th 2015
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
Threads 1
Posts 5
1:21:05 AM PDT - Thu, Apr 23rd 2015
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
Threads 1
Posts 5
5:49:33 AM PDT - Mon, Apr 27th 2015
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
Threads 3
Posts 855
3:46:00 PM PDT - Mon, Apr 27th 2015
It seems reasonable to me.
Does it seem to produce meaningful results?

Clicked A Few Times
Threads 1
Posts 5
7:02:27 AM PDT - Tue, Apr 28th 2015
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
Threads 3
Posts 855
2:16:45 PM PDT - Tue, Apr 28th 2015
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
Threads 1
Posts 5
5:38:30 AM PDT - Wed, Apr 29th 2015
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
Threads 3
Posts 855
9:27:19 AM PDT - Wed, Apr 29th 2015
Thanks for your contribution


Forum >> NWChem's corner >> General Topics



Who's here now Members 0 Guests 1 Bots/Crawler 0

AWC's: 2.5.10 MediaWiki - Stand Alone Forum Extension
Forum theme style by: AWC
Retrieved from "http://www.nwchem-sw.org/index.php/Special:AWCforum/st/id1659/Accessing_electron_density_from_....html"