Merge pull request #4490 from aannabe/density_fix

Fix units in qdens-radial and add documentation

docs/analyzing.rst

@@ -1352,3 +1352,105 @@ sufficient number of blocks remaining following any reblocking.
 
 When twist averaging, the group tag (e.g. ``g000`` or similar) will be
 replaced with ``avg`` in the names of the outputted files.
+
+.. _qdens-radial:
+
+Using the qdens-radial tool to estimate atomic occupations
+----------------------------------------------------------
+
+Once the ``*Density*.xsf`` files are obtained from ``qdens``, one can use the ``qdens-radial`` tool to calculate the on-site populations.
+Given a set of species and radii (in Angstroms), this tool will generate a radial density – angular average – around the atomic sites up to the specified radius.
+The radial density can be chosen to be non-cumulative or cumulative (integrated).
+
+::
+
+  >qdens-radial
+
+  Usage: qdens-radial [options] xsf_file
+  
+  Options:
+    --version             show program's version number and exit
+    -h, --help            Print help information and exit (default=False).
+    -v, --verbose         Print detailed information (default=False).
+    -r RADII, --radii=RADII
+                          List of cutoff radii (default=None).
+    -s SPECIES, --species=SPECIES
+                          List of species (default=None).
+    -a APP, --app=APP     Source that generated the .xsf file. Options: "pwscf",
+                          "qmcpack" (default=qmcpack).
+    -c, --cumulative      Evaluate cumulative radial density at cutoff radii
+                          (default=False).
+    --vmc=VMC_FILE        Location of VMC to be used for extrapolating mixed-
+                          estimator bias (default=None).
+    --write               Write extrapolated values to qmc-extrap.xsf
+                          (default=False).
+    --vmcerr=VMC_ERR_FILE
+                          Location of VMC+err to be used for extrapolating
+                          mixed-estimator bias and resampling (default=None).
+    --dmcerr=DMC_ERR_FILE
+                          Location of DMC+err to be used for resampling
+                          (default=None).
+    --seed=RANDOM_SEED    Random seed used for re-sampling. (default=None).
+    -n NSAMPLES, --nsamples=NSAMPLES
+                          Number of samples for resampling (default=50).
+    -p, --plot            Show plots interactively (default=False).
+
+
+Usage examples
+~~~~~~~~~~~~~~
+
+Below are example use cases for the H2O molecule using DMC data.
+
+Plot DMC non-cumulative radial density of the O atom:
+::
+  qdens-radial -p -s O -r 1 dmc.s002.Density_q.xsf
+
+Plot DMC cumulative radial density of the O atom:
+::
+  qdens-radial -p -s O -r 1 -c dmc.s002.Density_q.xsf
+
+For the cumulative case, ``qdens-radial`` will also print the cumulative value at the specified radius, i.e., an estimate of the atomic occupation.
+
+Estimate of the DMC atomic occupation:
+::
+  qdens-radial -p -s O -r 1.1 -c dmc.s002.Density_q.xsf
+
+Output:
+::
+  Cumulative Value of O Species at Cutoff 1.1 is: 6.55517033828574
+
+
+One can also get an extrapolated estimate (mixed-estimator bias) for this quantity by providing a VMC ``.xsf`` file.
+
+Estimate of the extrapolated atomic occupation:
+::
+  qdens-radial -p -s O -r 1.1 -c --vmc=dmc.s000.Density_q.xsf dmc.s002.Density_q.xsf
+  
+Output:
+::
+  Extrapolating from VMC and DMC densities...
+  Cumulative Value of O Species at Cutoff 1.1 is: 6.576918233167152
+
+Error bars
+~~~~~~~~~~
+
+One can "resample" the density at each grid point to obtain an estimate of the error bar. Recipe:
+
+1. Use error bars from ``*.Density_q+err.xsf`` file and draw samples from a Gaussian
+distribution with a standard deviation that matches the error bar.
+2. Calculate occupations with resampled data and calculate standard deviation
+to obtain the error bar on the occupation.
+3. Make sure the number of samples (``-n``) is converged.
+
+Estimate DMC atomic occupation with error bar:
+::
+  qdens-radial -p -s O -r 1.1 -c -n 20 --dmcerr=dmc.s002.Density_q+err.xsf dmc.s002.Density_q.xsf
+
+
+Output:
+::
+  Resampling to obtain error bar (NOTE: This can be slow)...
+  Will compute 20 samples...
+  ...
+  Cumulative Value of O Species at Cutoff 1.1 is: 6.55517033828574+/-0.001558553749396279
+

nexus/bin/qdens-radial

@@ -381,9 +381,11 @@ class QMCDensityRadialProcessor(QDRBase):
                 cds.read_xsf(resample_xsf_data)
                 equiv_atoms = list(cds.info.structure.equivalent_atoms().keys())
                 if opt.app=='pwscf':
-                    cds.change_distance_units('B')
+                    cds.set_attribute('density_units','B') # Quantum Espresso uses (e/Bohr^3) units for charge and spin density
+                    cds.change_density_units('A') # Change to (e/Angstrom^3) units for consistency with QMCPACK case
                 else:
                     cds.volume_normalize()
+                    cds.set_attribute('density_units','A') # xsf format uses Angstrom units which makes the QMCPACK density units (e/Angstrom^3) after volume normalization
                 #end if
 
                 if i==0:
@@ -456,9 +458,11 @@ class QMCDensityRadialProcessor(QDRBase):
         #end if
 
         if opt.app=='pwscf':
-            cd.change_distance_units('B')
+            cd.set_attribute('density_units','B') # Quantum Espresso uses (e/Bohr^3) units for charge and spin density
+            cd.change_density_units('A') # Change to (e/Angstrom^3) units for consistency with QMCPACK case
         else:
             cd.volume_normalize()
+            cd.set_attribute('density_units','A') # xsf format uses Angstrom units which makes the QMCPACK density units (e/Angstrom^3) after volume normalization
         #end if
 
         self.log('\nNorm:',cd.norm())

nexus/lib/observables.py

@@ -1191,22 +1191,20 @@ def norm(self,component=None):
 
     def change_distance_units(self,units):
         units_old = self.get_attribute('distance_units')
-        rscale    = 1.0/convert(1.0,units_old,units)
-        dscale    = 1./rscale**3
+        rscale    = convert(1.0,units_old,units)
         grid      = self.get_attribute('grid')
         grid.points *= rscale
-        for c in self.components():
-            c.values *= dscale
-        #end for
+        self.set_attribute('distance_units',units) # Update the object info to reflect the conversion
     #end def change_distance_units
 
 
     def change_density_units(self,units):
         units_old = self.get_attribute('density_units')
         dscale    = 1.0/convert(1.0,units_old,units)
         for c in self.components():
-            c.values *= dscale
+            c.values *= dscale**3
         #end for
+        self.set_attribute('density_units',units) # Update the object info to reflect the conversion
     #end def change_density_units
 
 
@@ -1339,6 +1337,7 @@ def plot_radial_density(self,component=None,show=True,cumulative=False,**kwargs)
         species = list(rdf.keys())
 
         dist_units = self.get_attribute('distance_units',None)
+        density_units = self.get_attribute('density_units',None)
 
         for cname in self.component_names:
             if cname in rdfs:
@@ -1353,10 +1352,14 @@ def plot_radial_density(self,component=None,show=True,cumulative=False,**kwargs)
                     #end if
                     plt.xlabel(xlabel)
                     if not cumulative:
-                        plt.ylabel('Radial density')
+                        ylabel = 'Radial density'
                     else:
-                        plt.ylabel('Cumulative radial density')
+                        ylabel = 'Cumulative radial density'
+                    #end if
+                    if density_units is not None:
+                        ylabel += ' (e/{}^3)'.format(density_units)
                     #end if
+                    plt.ylabel(ylabel)
                     plt.title('{} {} density'.format(s,cname))
                 #end for
             #end if

nexus/sphinx_docs/user-scripts.rst

@@ -530,7 +530,7 @@ finished.
 .. _user-limit-queue:
 
 Limiting the number of submitted jobs
--------------
+-------------------------------------
 Nexus will submit all eligible jobs at the same time unless told otherwise. This can be a large number when many calculations are present within the
 same project, e.g. various geometries or twists.  While this is fine on local resources, it might break the rules at computing centers such as ALCF
 where only 20 jobs can be submitted at the same time. In such cases, it is possible to specify the the size of the queue in Nexus to avoid monopolizing