Using the qmca tool to obtain total energies and related quantities¶
The qmca tool is the primary means of analyzing scalar-valued data
generated by QMCPACK. Output files that contain scalar-valued data are
*.scalar.dat and *.dmc.dat (see Output Overview for a detailed description of these files).
Quantities that are available for analysis in *.scalar.dat files
include the local energy and its variance, kinetic energy, potential
energy and its components, acceptance ratio, and the average CPU time
spent per block, among others. The *.dmc.dat files provide
information regarding the DMC walker population in addition to the local
energy.
Basic capabilities of qmca include calculating mean values and
associated error bars, processing multiple files at once in batched
fashion, performing twist averaging, plotting mean values by series, and
plotting traces (per block or step) of the underlying data. These
capabilities are explained with accompanying examples in the following
subsections.
To use qmca, installations of Python and NumPy must be present on
the local machine. For graphical plotting, the matplotlib module must
also be available.
An overview of all supported input flags to qmca can be obtained by
typing qmca at the command line with no other inputs (also try
qmca-x for a short list of examples):
>qmca
no files provided, please see help info below
Usage: qmca [options] [file(s)]
Options:
--version show program's version number and exit
-v, --verbose Print detailed information (default=False).
-q QUANTITIES, --quantities=QUANTITIES
Quantity or list of quantities to analyze. See names
and abbreviations below (default=all).
-u UNITS, --units=UNITS
Desired energy units. Can be Ha (Hartree), Ry
(Rydberg), eV (electron volts), kJ_mol (k.
joule/mole), K (Kelvin), J (Joules) (default=Ha).
-e EQUILIBRATION, --equilibration=EQUILIBRATION
Equilibration length in blocks (default=auto).
-a, --average Average over files in each series (default=False).
-w WEIGHTS, --weights=WEIGHTS
List of weights for averaging (default=None).
-b, --reblock (pending) Use reblocking to calculate statistics
(default=False).
-p, --plot Plot quantities vs. series (default=False).
-t, --trace Plot a trace of quantities (default=False).
-h, --histogram (pending) Plot a histogram of quantities
(default=False).
-o, --overlay Overlay plots (default=False).
--legend=LEGEND Placement of legend. None for no legend, outside for
outside legend (default=upper right).
--noautocorr Do not calculate autocorrelation. Warning: error bars
are no longer valid! (default=False).
--noac Alias for --noautocorr (default=False).
--sac Show autocorrelation of sample data (default=False).
--sv Show variance of sample data (default=False).
-i, --image (pending) Save image files (default=False).
-r, --report (pending) Write a report (default=False).
-s, --show_options Print user provided options (default=False).
-x, --examples Print examples and exit (default=False).
--help Print help information and exit (default=False).
-d DESIRED_ERROR, --desired_error=DESIRED_ERROR
Show number of samples needed for desired error bar
(default=none).
-n PARTICLE_NUMBER, --enlarge_system=PARTICLE_NUMBER
Show number of samples needed to maintain error bar on
larger system: desired particle number first, current
particle number second (default=none)
Obtaining a statistically correct mean and error bar¶
A rough guess at the mean and error bar of the local energy can be
obtained in the following way with qmca:
In this case the VMC energy of an 8-atom cell of diamond is estimated to
be \(-45.876(2)\) Hartrees (Ha). This rough guess should not be used
for production-level or publication-quality estimates.
To obtain production-level results, the underlying data should first be
inspected visually to ensure that all data included in the averaging can
be attributed to a distribution sharing the same mean. The first steps
of essentially any MC calculation (the “equilibration phase”) do not
belong to the equilibrium distribution and should be excluded from
estimates of the mean and its error bar.
We can plot a data trace (-t) of the local energy in the
following way:
>qmca-t-qe-e0qmc.s000.scalar.dat
The -e0 part indicates that we do not want any data
to be initially excluded from the calculation of averages. The resulting
plot is shown in Fig. 4. The unphysical
equilibration period is visible on the left side of the plot.
Most of the data fluctuates around a well-defined mean (consistent
variations around a flat line). This property is important to verify
by plotting the trace for each QMC run.
If we exclude none of the equilibration data points, we get an
erroneous estimate of \(-45.870(2)\) Ha for the local energy:
The equilibration period is typically estimated by eye, though a few conservative values should be checked to ensure that the mean remains
unaffected. In this dataset, the equilibration appears to have been
reached after 100 or so samples. After excluding the first 100
VMC blocks from the analysis we get
This estimate (\(-45.877(2)\) Ha) differs significantly from the
\(-45.870(2)\) Ha figure obtained from the full set of data, but it
agrees with the rough estimate of \(-45.876(2)\) Ha obtained with
the abbreviated command (qmca-qeqmc.s000.scalar.dat). This is
because qmca makes a heuristic guess at the equilibration period and
got it reasonably correct in this case. In many cases, the heuristic
guess fails and should not be relied on for quality results.
We have so far obtained a statistically correct mean. To obtain a
statistically correct error bar, it is best to include \(\sim\)100
or more statistically independent samples. An estimate of the number of
independent samples can be obtained by considering the autocorrelation
time, which is essentially a measure of the number of samples that must
be traversed before an uncorrelated/independent sample is reached. We
can get an estimate of the autocorrelation time in the following way:
The flag –sac stands for (s)how (a)uto(c)orrelation. In this case,
the autocorrelation estimate is \(4.8\approx 5\) samples. Since the
total run contained 800 samples and we have excluded 100 of them, we can
estimate the number of independent samples as \((800-100)/5=140\).
In this case, the error bar is expected to be estimated reasonably well.
Keep in mind that the error bar represents the expected range of the
mean with a certainty of only \(\sim 70\%\); i.e., it is a one sigma
error bar. The actual mean value will lie outside the range indicated by
the error bar in 1 out of every 3 runs, and in a set of 20 runs 1 value
can be expected to deviate from its estimate by twice the error bar.
Wavefunction optimization is a highly nonlinear and sometimes
sensitive process. As such, there is a risk that systematic
errors encountered at this stage of the QMC process can be propagated
into subsequent (expensive) DMC runs unless they are guarded against
with vigilance.
In this section we again consider an 8-atom cell of diamond but now in
the context of Jastrow optimization (one- and two-body terms). In
optimization runs it is often preferable to use a large number of
warmupsteps (\(\sim 100\)) so that equilibration bias does not
propagate into the optimization process. We can check that the added
warm-up has had its intended effect by again checking the local energy
trace:
>qmca-t-qe*scalar*
The resulting plot can be found in Fig. 5. In
this case sufficient warmupsteps were used to exit the equilibration
period before samples were collected and we can proceed without using
the -e option with qmca.
After inspecting the trace, we should inspect the text output from
qmca, now including the total energy and its variance:
The flags -qev requested the energy (e) and the variance
(v). For this combination of quantities, a third column (ratio)
is printed containing the ratio of the variance and the absolute value
of the local energy. The variance/energy ratio is an intensive quantity
and is useful to inspect regardless of the system under study.
Successful optimization of molecules and solids of any size generally
result in comparable values for the variance/energy ratio.
The first line of the output (series0) corresponds to the local
energy and variance of the system without a Jastrow factor (all Jastrow
coefficients were initialized to zero in this case), reflecting the
quality of the orbitals alone. For pseudopotential systems, a
variance/energy ratio \(>0.20\) Ha generally indicates there is a
problem with the input orbitals that needs to be resolved before
performing wavefunction optimization.
The subsequent lines correspond to energies and variances of
intermediate parameterizations of the trial wavefunction during the
optimization process. The output line containing optseries1, for
example, corresponds to the trial wavefunction parameterized during the
series0 step (the parameters of this wavefunction would be found in
an output file matching *s000*opt.xml). The first thing to check
about the resulting optimization is again the variance/energy ratio. For
pseudopotential systems, a variance/energy ratio \(<0.03\) Ha is
consistent with a trial wavefunction of production quality, and values
of \(0.01\) Ha are rarely obtainable for standard Slater-Jastrow
wavefunctions. By this metric, all parameterizations obtained for
optimizations performed in series 0-6 are of comparable quality (note
that the quality of the wavefunction obtained during optimization series
7 is effectively unknown).
A good way to further discriminate among the parameterizations is to
plot the energy and variance as a function of series with qmca:
>qmca-p-qevopt*scalar.dat
The -p option results in plots of means plus error bars
vs. series for all requested quantities.
The resulting plots for the local energy and variance are shown
in Fig. 6. In this case, the resulting energies
and variances are statistically indistinguishable for all optimization
cycles.
A good way to choose the optimal wavefunction for use in DMC is to
select the one with the lowest statistically significant energy within
the set of optimized wavefunctions with reasonable variance (e.g., among
those with a variance/energy ratio \(<0.03\) Ha). For
pseudopotential calculations, minimizing according to the total energy
is recommended to reduce locality errors in DMC.
Judging the quality of the DMC projection process requires more
care than is needed in VMC. To reduce bias, a small
time step is required in the approximate projector but this also
leads to slow equilibration and long autocorrelation times.
Systematic errors in the projection process can also arise from
statistical fluctuations due to pseudopotentials or from trial
wavefunctions with larger-than-necessary variance.
To illustrate the problems that can arise with respect to slow
equilibration and long autocorrelation times, we consider the 8-atom
diamond system with VMC (\(200\) blocks of \(160\) steps)
followed by DMC (\(400\) blocks of \(5\) steps) with a small
time step (\(0.002\) Ha\(^{-1}\)). A good first step in
assessing the quality of any DMC run is to plot the trace of the local
energy:
>qmca-t-qe-e0*scalar*
The resulting trace plot is shown in Fig. 7. As
always, the DMC local energy decreases exponentially away from the VMC
value, but in this case it takes a long time to do so. At least half of
the DMC run is inefficiently consumed by equilibration. If we are not
careful to inspect and remove the transient, the estimated DMC energy
will be strongly biased by the transient as shown by the horizontal red
line (estimated mean) in the figure. The autocorrelation time is also
large (\(\sim 12\) blocks):
Of the included 200 blocks, fewer than 20 contribute to the estimated error
bar, indicating that we cannot trust the reported error bar.
This can also be demonstrated directly from the data. If we halve the number
of included samples to 100, we expect from Gaussian statistics
that the error bar will grow by a factor of \(\sqrt{2}\), but instead we
get
which erroneously shows an estimated increase in the error bar by a
factor of about 2. Overall, this run is simply too short to gain
meaningful information.
Consider the case in which we are interested in the cohesive energy of
diamond, and, after having performed a time step study of the cohesive
energy, we have found that the energy difference between bulk diamond
and atomic carbon converges to our required accuracy with a larger time
step of \(0.01\) Ha\(^{-1}\). In a production setting, a small
cell could be used to determine the appropriate time step, while a
larger cell would subsequently be used to obtain a converged cohesive
energy, though for purposes of demonstration we still proceed here with
the 8-atom cell. The new time step of \(0.01\) Ha\(^{-1}\)
will result in a shorter autocorrelation time than the smaller time step
used previously, but we would like to shorten the equilibration time
further still. This can be achieved by using a larger time step (say
\(0.02\) Ha\(^{-1}\)) in a short intermediate DMC run used to
walk down the transient. The rapidly achieved equilibrium with the
\(0.02\) Ha\(^{-1}\) time step projector will be much nearer
to the \(0.01\) Ha\(^{-1}\) time step we seek than the
original VMC equilibrium, so we can expect a shortened secondary
equilibration time in the production \(0.01\) Ha\(^{-1}\) time
step run. Note that this procedure is fully general, even if having to
deal with an even shorter time step (e.g., \(0.002\)
Ha\(^{-1}\)) for a particular problem.
We now rerun the previous example but with an intermediate DMC
calculation using \(40\) blocks of \(5\) steps with a time step
of \(0.02\) Ha\(^{-1}\), followed by a production DMC
calculation using \(400\) blocks of \(10\) steps with a time
step of \(0.01\) Ha\(^{-1}\). We again plot the local energy
trace using qmca:
>qmca-t-qe-e0*scalar*
with the result shown in Fig. 8.
The projection transient has been effectively contained in the
short DMC run with a larger time step. As expected, the
production run contains only a short equilibration period.
Removing the first 20 blocks as a precaution, we obtain an estimate
of the total energy in VMC and DMC:
Notice that the variance/energy ratio in DMC (\(0.220\) Ha) is
similar to but slightly smaller than that obtained with VMC
(\(0.235\) Ha). If the DMC variance/energy ratio is ever
significantly larger than with VMC, this is cause to be concerned about
the correctness of the DMC run. Also notice the estimated
autocorrelation time (\(\sim 5\) blocks). This leaves us with an
estimated \(\sim 76\) independent samples, though we should recall
that the autocorrelation time is also a statistical estimate that can be
improved with more data. We can gain a better estimate of the
autocorrelation time by using the *.dmc.dat files, which contain
output data resolved per step rather than per block (there are
\(10\times\) more steps than blocks in this example case):
This results in an estimated autocorrelation time of \(\sim 31\)
steps, or \(\sim 3\) blocks, indicating that we actually have
\(\sim 122\) independent samples, which should be sufficient to
obtain a trustworthy error bar. Our final DMC total energy is estimated
to be \(-46.0329(2)\) Ha.
Another simulation property that should be explicitly monitored
is the behavior of the DMC walker population. Data regarding the
walker population is contained in the *.dmc.dat files.
In Fig. 9 we show the trace of the DMC
walker population for the current run:
Following a DMC run, the walker population should be checked for two
qualities: (1) that the population is sufficiently large (a number
\(>2,000\) is generally sufficient to reduce population control
bias) and (2) that the population fluctuates benignly around its
intended target value. In this case the target walker count (provided in
the input file) was \(2,048\) and we can confirm from the plot that
the population is simply fluctuating around this value. Also, from the
text output we have a dynamic population estimate of 2,050(5) walkers.
Rapid population reductions or increases—population explosions—are
indicative of problems with a run. These issues sometimes result from
using a considerably poor wavefunction (see comments regarding
variance/energy ratio in the preceding subsections). QMCPACK has
internal guards in place that prevent the population from exceeding
certain maximum and minimum bounds, so in particularly faulty runs one
might see the population “stabilize” to a constant value much larger or
smaller than the target. In such cases the cause(s) for the divergent
population behavior needs to be investigated and resolved before
proceeding further.
A number of other scalar-valued quantities are available with qmca.
To obtain text output for all quantities available, simply exclude the
-q option used in previous examples. The following example shows
output for a DMC calculation of the 8-atom diamond system from the
scalar.dat file:
See the output overview for scalar.dat
(The .scalar.dat file) and dmc.dat
(The .dmc.dat file) for more information about
these quantities. The data analysis aspects for these
quantities are essentially the same as for the local
energy as covered in the preceding subsections.
Quantities that do not belong to an equilibrium distribution
(e.g., BlockCPU) are somewhat different, though they
still exhibit statistical fluctuations.
Batch file processing is a common use case for qmca. If we consider
an “equation-of-state” calculation involving the 8-atom diamond cell we
have used so far, we might be interested in the total energy for the
various supercell volumes along the trajectory from compression to
expansion. After checking the traces
(qmca-t-qescale_*/vmc/*scalar*) to settle on a sensible
equilibration cutoff as discussed in the preceding subsections, we can
obtain the total energies all at once:
In this case, we are using a Jastrow factor optimized only at the
equilibrium geometry (scale_1.00) but with radial cutoffs restricted
to the Wigner-Seitz radius of the most compressed supercell
(scale_0.80) to avoid introducing wavefunction cusps at the cell
boundary (had we tried, QMCPACK would have aborted with a warning in
this case). It is clear that this restricted Jastrow factor is not an
optimal choice because it yields variance/energy ratios between
\(0.036\) and \(0.057\) Ha. This issue is largely a result of
our undersized (8-atom) supercell; larger cells should always be used in
real production calculations.
Batch processing is also possible for multiple quantities. If multiple
quantities are requested, an additional line is inserted to separate
results from different runs:
>qmca-q'e bc eff'-e40scale_*/vmc/*scalar*scale_0.80/vmc/qmcseries0LocalEnergy=-44.6710+/-0.0061BlockCPU=0.02986+/-0.00038Efficiency=38104.00+/-0.00scale_0.82/vmc/qmcseries0LocalEnergy=-44.9828+/-0.0058BlockCPU=0.02826+/-0.00013Efficiency=44483.91+/-0.00scale_0.84/vmc/qmcseries0LocalEnergy=-45.2283+/-0.0054BlockCPU=0.02747+/-0.00030Efficiency=52525.12+/-0.00scale_0.86/vmc/qmcseries0LocalEnergy=-45.4158+/-0.0055BlockCPU=0.02679+/-0.00013Efficiency=50811.55+/-0.00scale_0.88/vmc/qmcseries0LocalEnergy=-45.5702+/-0.0047BlockCPU=0.02598+/-0.00015Efficiency=74148.79+/-0.00scale_0.90/vmc/qmcseries0LocalEnergy=-45.6837+/-0.0050BlockCPU=0.02527+/-0.00011Efficiency=65714.98+/-0.00...
Twist averaging can be performed straightforwardly for any output
quantity listed in Obtaining other quantities with qmca.
We illustrate these capabilities by repeating the 8-atom diamond DMC
runs performed in Section Judging diffusion Monte Carlo runs at 8 real-valued
supercell twist angles (a \(2\times 2\times 2\) Monkhorst-Pack grid
centered at the \(\Gamma\)-point). Data traces for each twist can be
overlapped on the same plot:
>qmca-to-qe-e'30 20 30'*scalar*--legendoutside
The -o option requests the plots to be overlapped; otherwise,
8 separate plots would be generated. The
equilibration input -e'302030' cuts out from
the analyzed data the first 30 blocks for series 0 (VMC),
20 blocks for series 1 (intermediate DMC), and 30 blocks for
series 2 (production DMC). The resulting plot is shown in
Fig. 10.
Twist averaging is performed by providing the -a
option. If provided on its own, uniform weights are applied
to each twist angle. To obtain a trace plot with twist averaging
enforced, use a command similar to the following:
>qmca-a-t-qe-e'30 20 30'*scalar*
The resulting plot is shown in Fig. 11. As
can be seen from the trace plot, the chosen equilibration lengths are
appropriate, and we proceed to obtain the twist-averaged total energy
from the scalar.dat files
The degenerate twists grouped by set are \(\{0\}\),
\(\{1,2,4\}\), \(\{3,5,6\}\), and \(\{7\}\).
Alternatively, the run could have been performed at the four
unique (irreducible) twist angles only. We will emulate this situation by
analyzing data for twists 0, 1, 3, and 7 only. In a production setting
with irreducibly weighted twists, the run would be performed on these twists
alone; we reuse the uniform twist data for illustration purposes only.
We can use qmca to perform twist averaging with different
weights applied to each twist:
yielding a total energy value of \(-45.874(1)\) Ha, in agreement with the
uniform weighted twist average performed previously.
The decision of whether or not to perform irreducible weighted twist
averaging should be made on the basis of efficiency. The relative
efficiency of irreducible vs. uniform weighted twist averaging
depends on the irreducible weights and the ratio of the lengths of
the available sampling and equilibration periods. A formula for
the relative efficiency of these two cases is derived and discussed
in more detail in Appendix A: Derivation of twist averaging efficiency.
When working with many files or files with many entries, qmca might
take a long time to produce plots. The time delay is actually due to the
autocorrelation time estimate used to calculate error bars. The
calculation time for the autocorrelation scales as
\(\mathcal{O}(M^2)\), with \(M\) being the number of statistical
samples. If you are interested only in plotting traces and not in the
estimated error bars, the autocorrelation time estimation can be turned
off with the –noac option:
>qmca-t-qe-e20--noacqmc.s002.scalar.dat
Note that the resulting error bars printed to the console will be
underestimated and are not meaningful. Do not use –noac in
conjunction with the -p plotting option as these plots are of no use
without meaningful error bars.
Inspect the trace plots (-t option) for any oddities in the data.
Typical behavior is a short equilibration period followed by benign
fluctuations around a clear mean value. There should not be any large
spikes in the data. This applies to all runs (VMC, optimization,
DMC, etc.).
Remove all equilibration steps (-e option) from the data by
inspecting the trace plot.
Check the quality of the orbitals (standalone Jastrow-less VMC or
sometimes the first scalar file produced during optimization) by
inspecting the variance/energy ratio qmca-qev*scalar*. For
pseudopotential systems without a Jastrow, the variance/energy ratio
should not exceed \(0.2\) Ha; otherwise, there is a problem with
the orbitals.
Check the quality of the optimized Jastrow factor by inspecting the
variance/energy ratio. For pseudopotential systems with a Jastrow,
the variance/energy ratio should not exceed \(0.04\) Ha for
pseudopotential systems. A good Jastrow is indicated by a
variance/energy ratio in the range of \(0.01-0.03\) Ha. A value
less than \(0.01\) Ha is difficult to achieve.
Confirm that the optimization has converged by plotting the energy
and variance vs. optimization series (qmca-p-qev*scalar*). Do
not assume that optimization has converged in only a few cycles. Use
at least 10 cycles with about 100,000 samples unless you already have
experience with the system in question.
Optimize Jastrow factors according to energy minimization to reduce
locality errors arising from the use of nonlocal pseudopotentials in
DMC. A good approach is to optimize with a few cycles of variance
minimization followed by several cycles of energy minimization.
Occasionally try optimizing with more samples and/or cycles to see if
improved results are obtained.
If using a B-spline representation of the orbitals, converge the VMC
energy and variance with respect to the mesh size (controlled via
meshfactor). This is best done in the presence of any Jastrow factor
to reduce noise. Consider using the hybrid LMTO representation of the
orbitals as this can reduce both the VMC/DMC variance and the DMC
time step error, in addition to saving memory.
Check the variance/energy ratio of all production VMC and DMC
calculations. In all cases, the DMC ratio should be slightly less
than the VMC ratio and both should abide the preceding guidelines,
i.e., the ratio should be less than \(0.04\) Ha for
pseudopotential systems. The production ratio should also be
consistent with what is observed during wavefunction optimization.
Be aware of population control bias in DMC. Run with a population of
\(\sim 2,000\) or greater. Occasionally repeat a run using a
larger population to explicitly confirm that population control bias
is small.
Check the stability of the DMC walker population by plotting the
trace of the population size (qmca-t-qnw*dmc.dat). Verify
that the average walker population is consistent with the requested
value provided in the input.
In DMC, perform a time step study to obtain either (1) extrapolated
results or (2) a time step for future production where an energy
difference shows convergence (e.g., a band gap or defect formation
energy). For pseudopotential systems, converged time steps for many
systems are in the range of \(0.002-0.01\) Ha\(^{-1}\), but
the actual converged time step must be explicitly checked.
In periodic systems, converge the total energy with respect to the
size of the twist/k-point grid. Results for smaller systems can
easily be transferred to larger ones (e.g., a
\(2 \times 2 \times 2\) twist grid in a
\(2 \times 2 \times 2\) tiled cell is equivalent to a
\(1 \times 1 \times 1\) twist grid in a
\(4 \times 4 \times 4\) tiled cell).
In periodic systems, perform finite-size extrapolation including two
body corrections (needed for cohesive energy/phase stability studies)
unless it can be shown that finite-size effects cancel for the energy
difference in question (e.g., some defect formation energies).
Using the qmc-fit tool for statistical time step extrapolation and curve fitting¶
The qmc-fit tool is used to provide statistical estimates of
curve-fitting parameters based on QMCPACK data. Although qmc-fit
will eventually support many types of fitted curves (e.g., Morse
potential binding curves and various equation-of-state fitting curves),
it is currently limited to estimating fitting parameters related to time
step extrapolation.
The qmc-fit tool obtains estimates of fitting parameter means and
associated error bars via the “jack-knife” technique. This technique is
a powerful and general tool to obtain meaningful error bars for any
quantity that is related in a nonlinear fashion to an underlying set of
statistical data. For this reason, we give a brief overview of the
jackknife technique before proceeding with usage instructions for the
qmc-fit tool.
Consider \(N\) statistical variables \(\{x_n\}_{n=1}^N\) that
have been outputted by one or more simulation runs. If we have \(M\)
samples of each of the \(N\) variables, then the mean values of each
these variables can be estimated in the standard way, that is,
\(\bar{x}_n\approx \tfrac{1}{M}\sum_{m=1}^Mx_{nm}\).
Suppose we are interested in \(P\) statistical quantities
\(\{y_p\}_{p=1}^P\) that are related to the original \(N\)
variables by a known multidimensional function \(F\):
The relationship implied by \(F\) is completely general. For
example, the \(\{x_n\}\) might be elements of a matrix with
\(\{y_p\}\) being the eigenvalues, or \(F\) might be a fitting
procedure for \(N\) energies at different time steps with \(P\)
fitting parameters. An approximate guess at the mean value of
\(\vec{y}\) can be obtained by evaluating \(F\) at the mean
value of \(\vec{x}\) (i.e. \(F(\bar{x}_1\ldots\bar{x}_N)\)), but
with this approach we have no way to estimate the statistical error bar
of any \(\bar{y}_p\).
In the jackknife procedure, the statistical variability intrinsic to the
underlying data \(\{x_n\}\) is used to obtain estimates of the mean
and error bar of \(\{y_p\}\). We first construct a new set of
\(x\) statistical data by taking the average over all samples but
one:
In this section, we use a 32-atom supercell of MnO as an example system
for time step extrapolation. Data for this system has been collected in
DMC using the following sequence of time steps:
\(0.04,~0.02,~0.01,~0.005,~0.0025,~0.00125\) Ha\(^{-1}\). For
a typical production pseudopotential study, time steps in the range of
\(0.02-0.002\) Ha\(^{-1}\) are usually sufficient and it is
recommended to increase the number of steps/blocks by a factor of two
when the time step is halved. To perform accurate statistical fitting,
we must first understand the equilibration and autocorrelation
properties of the inputted local energy data. After plotting the local
energy traces (qmca-t-qe-e0./qmc*/*scalar*), it is clear that
an equilibration period of \(30\) blocks is reasonable. Approximate
autocorrelation lengths are also obtained with qmca:
The input arguments are as follows: ts indicates we are performing a
time step fit, -e30 is the equilibration period removed from each
set of scalar data, -b4 indicates the data will be reblocked by a
factor of 4 (e.g., a file containing 400 entries will be block averaged
into a new set of 100 before jackknife fitting), -s2 indicates that
the time step data begins with series 2 (scalar files matching
*s000* or *s001* are to be excluded), and -t ‘0.00125 0.0025
0.005 0.01 0.02 0.04’ provides a list of time step values corresponding
to the inputted scalar files. The -e and -b options can receive
a list of file-specific values (same format as -t) if desired. As
can be seen from the text output, the parameters for the linear fit are
printed with error bars obtained with jackknife resampling and the zero
time step “intercept” is \(-3848.19(4)\) Ha. In addition to text
output, the previous command will result in a plot of the fit with the
zero time step value shown as a red dot, as shown in the top panel of
Fig. 12.
Different fitting functions are supported via the -f option.
Currently supported options include linear (\(a+bt\)),
quadratic (\(a+bt+ct^2\)), and sqrt
(\(a+b\sqrt{t}+ct\)). Results for a quadratic fit are shown
subsequently and in the bottom panel of Fig. 12.
In this case, we find a zero time step estimate of \(-3848.25(5)\)
Ha\(^{-1}\). A time step of \(0.04\) Ha\(^{-1}\) might
be on the large side to include in time step extrapolation, and it is
likely to have an outsize influence in the case of linear extrapolation.
Upon excluding this point, linear extrapolation yields a zero timestep
value of \(-3848.22(4)\) Ha\(^{-1}\). Note that quadratic
extrapolation can result in intrinsically larger uncertainty in the
extrapolated value. For example, when the \(0.04\) Ha\(^{-1}\)
point is excluded, the uncertainty grows by 50% and we obtain an
estimated value of \(-3848.28(7)\) instead.
Using the qdens tool to obtain electron densities¶
The qdens tool is provided to post-process the heavy density data
produced by QMCPACK and output the mean density (with and without
errorbars) in file formats viewable with, e.g., XCrysDen or VESTA. The
tool currently works only with the SpinDensity estimator in QMCPACK.
Note: this tool is provisional and may be changed or replaced at any
time. The planned successor to this tool (qstat) will expand access
to other observables and will retain at least the non-plotting
capabilities of qdens.
To use qdens, Nexus must be installed along with NumPy and H5Py. A
short list of example use cases are covered in the next section. Current
input flags are:
>qdensUsage:qdens[options][file(s)]Options:--versionshowprogram's version number and exit-h,--helpPrinthelpinformationandexit(default=False).-v,--verbosePrintdetailedinformation(default=False).-fFORMATS,--formats=FORMATSFormatorlistofformatsfordensityfileoutput.Options:dat,xsf,chgcar(default=None).-eEQUILIBRATION,--equilibration=EQUILIBRATIONEquilibrationlengthinblocks(default=0).-rREBLOCK,--reblock=REBLOCKBlockcoarseningfactor;useestimatedautocorrelationlength(default=None).-a,--averageAverageoverfilesineachseries(default=False).-wWEIGHTS,--weights=WEIGHTSListofweightsforaveraging(default=None).-iINPUT,--input=INPUTQMCPACKinputfilecontainingstructureandgridinformation(default=None).-sSTRUCTURE,--structure=STRUCTUREFilecontainingatomicstructure(default=None).-gGRID,--grid=GRIDDensitygriddimensions(default=None).-cCELL,--cell=CELLSimulationcellaxes(default=None).--lineplot=LINEPLOTProducealineplotalongtheselecteddimension:0,1,or2(default=None).--noplotDonotshowplotsinteractively(default=False).
Process a single file, excluding the first 40 blocks, and produce XSF files:
qdens-v-e40-fxsf-iqmc.in.xmlqmc.s000.stat.h5
Process files for all available series:
qdens-v-e40-fxsf-iqmc.in.xml*stat.h5
Combine groups of 10 adjacent statistical blocks together (appropriate if the
estimated autocorrelation time is about 10 blocks):
qdens-v-e40-r10-fxsf-iqmc.in.xmlqmc.s000.stat.h5
Apply different equilibration lengths and reblocking factors to each
series (below is appropriate if there are three series, e.g. s000,
s001, and s002):
Files postfixed with u relate to the up electron density, d to
down, u+d to the total charge density, and u-d to the difference
between up and down electron densities.
Files without err in the name contain only the mean, whereas files
with +err/-err in the name contain the mean plus/minus the
estimated error bar. Please use caution in interpreting the error bars
as their accuracy depends crucially on a correct estimation of the
autocorrelation time by the user (see -r option) and having a
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.