From a39c7a8795739dac61799a5af13cd07b38e7d2d3 Mon Sep 17 00:00:00 2001 From: Gani Date: Wed, 1 Mar 2023 15:28:28 -0500 Subject: [PATCH 1/3] qdens-radial fix units --- nexus/bin/qdens-radial | 8 ++++++-- nexus/lib/observables.py | 19 +++++++++++-------- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/nexus/bin/qdens-radial b/nexus/bin/qdens-radial index 2ccadec631..dcf7c265a5 100755 --- a/nexus/bin/qdens-radial +++ b/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()) diff --git a/nexus/lib/observables.py b/nexus/lib/observables.py index 5a32406e92..206407abc0 100644 --- a/nexus/lib/observables.py +++ b/nexus/lib/observables.py @@ -1191,13 +1191,10 @@ 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 refect the conversion #end def change_distance_units @@ -1205,8 +1202,9 @@ 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 refect 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 From 89fffbf05b10cd4907722b8d4af29d169152d8b6 Mon Sep 17 00:00:00 2001 From: Gani Date: Wed, 1 Mar 2023 17:12:28 -0500 Subject: [PATCH 2/3] documentation and fix typo --- docs/analyzing.rst | 102 +++++++++++++++++++++++++++++++++++++++ nexus/lib/observables.py | 4 +- 2 files changed, 104 insertions(+), 2 deletions(-) diff --git a/docs/analyzing.rst b/docs/analyzing.rst index ce241c8246..91b5069dc9 100644 --- a/docs/analyzing.rst +++ b/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 + diff --git a/nexus/lib/observables.py b/nexus/lib/observables.py index 206407abc0..b1eb55e031 100644 --- a/nexus/lib/observables.py +++ b/nexus/lib/observables.py @@ -1194,7 +1194,7 @@ def change_distance_units(self,units): rscale = convert(1.0,units_old,units) grid = self.get_attribute('grid') grid.points *= rscale - self.set_attribute('distance_units',units) # Update the object info to refect the conversion + self.set_attribute('distance_units',units) # Update the object info to reflect the conversion #end def change_distance_units @@ -1204,7 +1204,7 @@ def change_density_units(self,units): for c in self.components(): c.values *= dscale**3 #end for - self.set_attribute('density_units',units) # Update the object info to refect the conversion + self.set_attribute('density_units',units) # Update the object info to reflect the conversion #end def change_density_units From 606c24af7c5b91e577ee098ee77bbd22680e9992 Mon Sep 17 00:00:00 2001 From: Ye Luo Date: Wed, 1 Mar 2023 19:44:50 -0600 Subject: [PATCH 3/3] Fix sphinx_docs warning. --- nexus/sphinx_docs/user-scripts.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nexus/sphinx_docs/user-scripts.rst b/nexus/sphinx_docs/user-scripts.rst index 7ede4e7e0d..c26e7de763 100644 --- a/nexus/sphinx_docs/user-scripts.rst +++ b/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