scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
qmcpack/docs/analyzing.rst

1603 lines
68 KiB
ReStructuredText
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

.. _analyzing:
Analyzing QMCPACK data
======================
.. _qmca:
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 :ref:`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)
.. _qmca-mean-error:
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``:
::
>qmca -q e qmc.s000.scalar.dat
qmc series 0 LocalEnergy = -45.876150 +/- 0.017688
In this case the VMC energy of an 8-atom cell of diamond is estimated to
be :math:`-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 -q e -e 0 qmc.s000.scalar.dat
The ``-e 0`` part indicates that we do not want any data
to be initially excluded from the calculation of averages. The resulting
plot is shown in :numref:`fig4`. The unphysical
equilibration period is visible on the left side of the plot.
.. _fig4:
.. figure:: /figs/qmca_mean_error_trace.png
:width: 400
:align: center
Trace of the VMC local energy for an 8-atom cell of diamond generated
with ``qmca``. The x-axis (“samples”) refers to the VMC block index in
this case.
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 :math:`-45.870(2)` Ha for the local energy:
::
>qmca -q e -e 0 qmc.s000.scalar.dat
qmc series 0 LocalEnergy = -45.870071 +/- 0.018072
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
::
>qmca -q e -e 100 qmc.s000.scalar.dat
qmc series 0 LocalEnergy = -45.877363 +/- 0.017432
This estimate (:math:`-45.877(2)` Ha) differs significantly from the
:math:`-45.870(2)` Ha figure obtained from the full set of data, but it
agrees with the rough estimate of :math:`-45.876(2)` Ha obtained with
the abbreviated command (``qmca -q e qmc.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 :math:`\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:
::
>qmca -q e -e 100 qmc.s000.scalar.dat --sac
qmc series 0 LocalEnergy = -45.877363 +/- 0.017432 4.8
The flag ``–sac`` stands for (s)how (a)uto(c)orrelation. In this case,
the autocorrelation estimate is :math:`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 :math:`(800-100)/5=140`.
In this case, the error bar is expected to be estimated reasonably well.
.. _fig5:
.. figure:: /figs/qmca_judge_opt.png
:width: 400
:align: center
Trace of the local energy during one- and two-body Jastrow optimizations
for an 8-atom cell of diamond generated with ``qmca``. Data for each
optimization cycle (QMCPACK series) is separated by a vertical black
line.
Keep in mind that the error bar represents the expected range of the
mean with a certainty of only :math:`\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.
.. qmca-judge-opt:
Judging wavefunction optimization
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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`` (:math:`\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 -q e *scalar*
The resulting plot can be found in :numref:`fig5`. 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:
::
>qmca -q ev opt*scalar.dat
LocalEnergy Variance ratio
opt series 0 -44.823616 +/- 0.007430 7.054219 +/- 0.041998 0.1574
opt series 1 -45.877643 +/- 0.003329 1.095362 +/- 0.041154 0.0239
opt series 2 -45.883191 +/- 0.004149 1.077942 +/- 0.021555 0.0235
opt series 3 -45.877524 +/- 0.003094 1.074047 +/- 0.010491 0.0234
opt series 4 -45.886062 +/- 0.003750 1.061707 +/- 0.014459 0.0231
opt series 5 -45.877668 +/- 0.003475 1.091585 +/- 0.021637 0.0238
opt series 6 -45.877109 +/- 0.003586 1.069205 +/- 0.009387 0.0233
opt series 7 -45.882563 +/- 0.004324 1.058771 +/- 0.008651 0.0231
The flags ``-q ev`` 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 (``series 0``) 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 :math:`>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 ``opt series 1``, for
example, corresponds to the trial wavefunction parameterized during the
``series 0`` 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 :math:`<0.03` Ha is
consistent with a trial wavefunction of production quality, and values
of :math:`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 -q ev opt*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 :numref:`fig6`. 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 :math:`<0.03` Ha). For
pseudopotential calculations, minimizing according to the total energy
is recommended to reduce locality errors in DMC.
.. image:: /figs/qmca_opt_energy.png
:width: 400
:align: center
.. _fig6:
.. figure:: /figs/qmca_opt_variance.png
:width: 400
:align: center
Energy and variance vs. optimization series for an 8-atom cell of
diamond as plotted by ``qmca``.
.. _qmca-judge-dmc:
Judging diffusion Monte Carlo runs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 (:math:`200` blocks of :math:`160` steps)
followed by DMC (:math:`400` blocks of :math:`5` steps) with a small
time step (:math:`0.002` Ha\ :math:`^{-1}`). A good first step in
assessing the quality of any DMC run is to plot the trace of the local
energy:
::
>qmca -t -q e -e 0 *scalar*
.. _fig7:
.. figure:: /figs/qmca_short_dmc.png
:width: 400
:align: center
Trace of the local energy for VMC followed by DMC with a small time step
(:math:`0.002` Ha\ :math:`^{-1}`) for an 8-atom cell of diamond
generated with ``qmca``.
The resulting trace plot is shown in :numref:`fig7`. 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 (:math:`\sim 12` blocks):
::
>qmca -q e -e 200 --sac *s001.scalar*
qmc series 1 LocalEnergy = -46.045720 +/- 0.004813 11.6
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 :math:`\sqrt{2}`, but instead we
get
::
>qmca -q e -e 300 *s001.scalar*
qmc series 1 LocalEnergy = -46.048537 +/- 0.009280
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 :math:`0.01` Ha\ :math:`^{-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 :math:`0.01` Ha\ :math:`^{-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
:math:`0.02` Ha\ :math:`^{-1}`) in a short intermediate DMC run used to
walk down the transient. The rapidly achieved equilibrium with the
:math:`0.02` Ha\ :math:`^{-1}` time step projector will be much nearer
to the :math:`0.01` Ha\ :math:`^{-1}` time step we seek than the
original VMC equilibrium, so we can expect a shortened secondary
equilibration time in the production :math:`0.01` Ha\ :math:`^{-1}` time
step run. Note that this procedure is fully general, even if having to
deal with an even shorter time step (e.g., :math:`0.002`
Ha\ :math:`^{-1}`) for a particular problem.
We now rerun the previous example but with an intermediate DMC
calculation using :math:`40` blocks of :math:`5` steps with a time step
of :math:`0.02` Ha\ :math:`^{-1}`, followed by a production DMC
calculation using :math:`400` blocks of :math:`10` steps with a time
step of :math:`0.01` Ha\ :math:`^{-1}`. We again plot the local energy
trace using ``qmca``:
::
>qmca -t -q e -e 0 *scalar*
with the result shown in :numref:`fig8`.
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:
::
>qmca -q ev -e 20 --sac qmc.*.scalar.dat
LocalEnergy Variance ratio
qmc series 0 -45.881042 +/- 0.001283 1.0 1.076726 +/- 0.007013 1.0 0.0235
qmc series 1 -46.040814 +/- 0.005046 3.9 1.011303 +/- 0.016807 1.1 0.0220
qmc series 2 -46.032960 +/- 0.002077 5.2 1.014940 +/- 0.002547 1.0 0.0220
.. _fig8:
.. figure:: /figs/qmca_accel_dmc.png
:width: 400
:align: center
Trace of the local energy for VMC followed by a short intermediate DMC
with a large time step (:math:`0.02` Ha\ :math:`^{-1}`) and finally a
production DMC run with a time step of :math:`0.01` Ha\ :math:`^{-1}`.
Calculations were performed in an 8-atom cell of diamond.
Notice that the variance/energy ratio in DMC (:math:`0.220` Ha) is
similar to but slightly smaller than that obtained with VMC
(:math:`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 (:math:`\sim 5` blocks). This leaves us with an
estimated :math:`\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
:math:`10\times` more steps than blocks in this example case):
::
>qmca -q ev -e 200 --sac qmc.s002.dmc.dat
LocalEnergy Variance ratio
qmc series 2 -46.032909 +/- 0.002068 31.2 1.015781 +/- 0.002536 1.4 0.0221
This results in an estimated autocorrelation time of :math:`\sim 31`
steps, or :math:`\sim 3` blocks, indicating that we actually have
:math:`\sim 122` independent samples, which should be sufficient to
obtain a trustworthy error bar. Our final DMC total energy is estimated
to be :math:`-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 :numref:`fig9` we show the trace of the DMC
walker population for the current run:
::
>qmca -t -q nw *dmc.dat
qmc series 1 NumOfWalkers = 2056.905405 +/- 8.775527
qmc series 2 NumOfWalkers = 2050.164160 +/- 4.954850
Following a DMC run, the walker population should be checked for two
qualities: (1) that the population is sufficiently large (a number
:math:`>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 :math:`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.
.. _fig9:
.. figure:: /figs/qmca_pop_trace.png
:width: 400
:align: center
Trace of the DMC walker population for an 8-atom cell of diamond
obtained with ``qmca``.
.. _qmca-other-quantities:
Obtaining other quantities
~~~~~~~~~~~~~~~~~~~~~~~~~~
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:
::
>qmca -e 20 qmc.s002.scalar.dat
qmc series 2
LocalEnergy = -46.0330 +/- 0.0021
Variance = 1.0149 +/- 0.0025
Kinetic = 33.851 +/- 0.019
LocalPotential = -79.884 +/- 0.020
ElecElec = -11.4483 +/- 0.0083
LocalECP = -22.615 +/- 0.029
NonLocalECP = 5.2815 +/- 0.0079
IonIon = -51.10 +/- 0.00
LocalEnergy_sq = 2120.05 +/- 0.19
BlockWeight = 20514.27 +/- 48.38
BlockCPU = 1.4890 +/- 0.0038
AcceptRatio = 0.9963954 +/- 0.0000055
Efficiency = 71.88 +/- 0.00
TotalTime = 565.80 +/- 0.00
TotalSamples = 7795421 +/- 0
Similarly, for the ``dmc.dat`` file we get
::
>qmca -e 20 qmc.s002.dmc.dat
qmc series 2
LocalEnergy = -46.0329 +/- 0.0020
Variance = 1.0162 +/- 0.0025
TotalSamples = 8201275 +/- 0
TrialEnergy = -46.0343 +/- 0.0023
DiffEff = 0.9939150 +/- 0.0000088
Weight = 2050.23 +/- 4.82
NumOfWalkers = 2050 +/- 5
LivingFraction = 0.996427 +/- 0.000021
AvgSentWalkers = 0.2625 +/- 0.0011
Any subset of desired quantities can be obtained by using the ``-q``
option with either the full names of the quantities just listed
::
>qmca -q 'LocalEnergy Kinetic LocalPotential' -e 20 qmc.s002.scalar.dat
qmc series 2
LocalEnergy = -46.0330 +/- 0.0021
Kinetic = 33.851 +/- 0.019
LocalPotential = -79.884 +/- 0.020
or with their corresponding abbreviations.
::
>qmca -q ekp -e 20 qmc.s002.scalar.dat
qmc series 2
LocalEnergy = -46.0330 +/- 0.0021
Kinetic = 33.851 +/- 0.019
LocalPotential = -79.884 +/- 0.020
Abbreviations for each quantity can be found by typing ``qmca`` at the
command line with no other input. This following is a current list:
::
Abbreviations and full names for quantities:
ar = AcceptRatio
bc = BlockCPU
bw = BlockWeight
ce = CorrectedEnergy
de = DiffEff
e = LocalEnergy
ee = ElecElec
eff = Efficiency
ii = IonIon
k = Kinetic
kc = KEcorr
l = LocalECP
le2 = LocalEnergy_sq
mpc = MPC
n = NonLocalECP
nw = NumOfWalkers
p = LocalPotential
sw = AvgSentWalkers
te = TrialEnergy
ts = TotalSamples
tt = TotalTime
v = Variance
w = Weight
See the output overview for ``scalar.dat``
(:ref:`scalardat-file`) and ``dmc.dat``
(:ref:`dmc-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.
.. _qmca-multiple-files:
Processing multiple files
~~~~~~~~~~~~~~~~~~~~~~~~~
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 -q e scale_*/vmc/*scalar*``) to settle on a sensible
equilibration cutoff as discussed in the preceding subsections, we can
obtain the total energies all at once:
.. code-block:: python
>qmca -q ev -e 40 scale_*/vmc/*scalar*
LocalEnergy Variance ratio
scale_0.80/vmc/qmc series 0 -44.670984 +/- 0.006051 2.542384 +/- 0.019902 0.0569
scale_0.82/vmc/qmc series 0 -44.982818 +/- 0.005757 2.413011 +/- 0.022626 0.0536
scale_0.84/vmc/qmc series 0 -45.228257 +/- 0.005374 2.258577 +/- 0.019322 0.0499
scale_0.86/vmc/qmc series 0 -45.415842 +/- 0.005532 2.204980 +/- 0.052978 0.0486
scale_0.88/vmc/qmc series 0 -45.570215 +/- 0.004651 2.061374 +/- 0.014359 0.0452
scale_0.90/vmc/qmc series 0 -45.683684 +/- 0.005009 1.988539 +/- 0.018267 0.0435
scale_0.92/vmc/qmc series 0 -45.751359 +/- 0.004928 1.913282 +/- 0.013998 0.0418
scale_0.94/vmc/qmc series 0 -45.791622 +/- 0.005026 1.843704 +/- 0.014460 0.0403
scale_0.96/vmc/qmc series 0 -45.809256 +/- 0.005053 1.829103 +/- 0.014536 0.0399
scale_0.98/vmc/qmc series 0 -45.806235 +/- 0.004963 1.775391 +/- 0.015199 0.0388
scale_1.00/vmc/qmc series 0 -45.783481 +/- 0.005293 1.726869 +/- 0.012001 0.0377
scale_1.02/vmc/qmc series 0 -45.741655 +/- 0.005627 1.681776 +/- 0.011496 0.0368
scale_1.04/vmc/qmc series 0 -45.685101 +/- 0.005353 1.682608 +/- 0.015423 0.0368
scale_1.06/vmc/qmc series 0 -45.615164 +/- 0.005978 1.652155 +/- 0.010945 0.0362
scale_1.08/vmc/qmc series 0 -45.543037 +/- 0.005191 1.646375 +/- 0.013446 0.0361
scale_1.10/vmc/qmc series 0 -45.450976 +/- 0.004794 1.707649 +/- 0.048186 0.0376
scale_1.12/vmc/qmc series 0 -45.371851 +/- 0.005103 1.686997 +/- 0.035920 0.0372
scale_1.14/vmc/qmc series 0 -45.265490 +/- 0.005311 1.631614 +/- 0.012381 0.0360
scale_1.16/vmc/qmc series 0 -45.161961 +/- 0.004868 1.656586 +/- 0.014788 0.0367
scale_1.18/vmc/qmc series 0 -45.062579 +/- 0.005971 1.671998 +/- 0.019942 0.0371
scale_1.20/vmc/qmc series 0 -44.960477 +/- 0.004888 1.651864 +/- 0.009756 0.0367
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
:math:`0.036` and :math:`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' -e 40 scale_*/vmc/*scalar*
scale_0.80/vmc/qmc series 0
LocalEnergy = -44.6710 +/- 0.0061
BlockCPU = 0.02986 +/- 0.00038
Efficiency = 38104.00 +/- 0.00
scale_0.82/vmc/qmc series 0
LocalEnergy = -44.9828 +/- 0.0058
BlockCPU = 0.02826 +/- 0.00013
Efficiency = 44483.91 +/- 0.00
scale_0.84/vmc/qmc series 0
LocalEnergy = -45.2283 +/- 0.0054
BlockCPU = 0.02747 +/- 0.00030
Efficiency = 52525.12 +/- 0.00
scale_0.86/vmc/qmc series 0
LocalEnergy = -45.4158 +/- 0.0055
BlockCPU = 0.02679 +/- 0.00013
Efficiency = 50811.55 +/- 0.00
scale_0.88/vmc/qmc series 0
LocalEnergy = -45.5702 +/- 0.0047
BlockCPU = 0.02598 +/- 0.00015
Efficiency = 74148.79 +/- 0.00
scale_0.90/vmc/qmc series 0
LocalEnergy = -45.6837 +/- 0.0050
BlockCPU = 0.02527 +/- 0.00011
Efficiency = 65714.98 +/- 0.00
...
.. _qmca-twist-average:
Twist averaging
~~~~~~~~~~~~~~~
Twist averaging can be performed straightforwardly for any output
quantity listed in :ref:`qmca-other-quantities` with ``qmca``.
We illustrate these capabilities by repeating the 8-atom diamond DMC
runs performed in Section :ref:`qmca-judge-dmc` at 8 real-valued
supercell twist angles (a :math:`2\times 2\times 2` Monkhorst-Pack grid
centered at the :math:`\Gamma`-point). Data traces for each twist can be
overlapped on the same plot:
::
>qmca -to -q e -e '30 20 30' *scalar* --legend outside
The ``-o`` option requests the plots to be overlapped; otherwise,
8 separate plots would be generated. The
equilibration input ``-e '30 20 30'`` 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
:numref:`fig10`.
.. _fig10:
.. figure:: /figs/qmca_twist_trace_overlap.png
:width: 400
:align: center
Overlapped energy traces from VMC to DMC for an 8-supercell diamond
obtained with ``qmca``. Data for each twist appears in a different
color.
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 -q e -e '30 20 30' *scalar*
The resulting plot is shown in :numref:`fig11`. 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
::
>qmca -a -q ev -e 30 --sac *s002.scalar*
LocalEnergy Variance ratio
avg series 2 -45.873369 +/- 0.000753 5.3 1.028751 +/- 0.001056 1.3 0.0224
and also from the ``dmc.dat`` files
::
>qmca -a -q ev -e 300 --sac *s002.dmc*
LocalEnergy Variance ratio
avg series 2 -45.873371 +/- 0.000741 30.5 1.028843 +/- 0.000972 1.6 0.0224
yielding a twist-averaged total energy of :math:`-45.8733(8)` Ha.
.. _fig11:
.. figure:: /figs/qmca_twist_average_trace.png
:width: 400
:align: center
Twist-averaged energy trace from VMC to DMC for an 8-supercell diamond obtained with ``qmca``.
As can be seen from :numref:`fig10`, some of the twist
angles are degenerate. This is seen more clearly in the text output
::
>qmca -q ev -e 30 *s002.scalar*
LocalEnergy Variance ratio
qmc.g000 series 2 -45.264510 +/- 0.001942 1.057065 +/- 0.002318 0.0234
qmc.g001 series 2 -46.035511 +/- 0.001806 1.015992 +/- 0.002836 0.0221
qmc.g002 series 2 -46.035410 +/- 0.001538 1.015039 +/- 0.002661 0.0220
qmc.g003 series 2 -46.047285 +/- 0.001898 1.018219 +/- 0.002588 0.0221
qmc.g004 series 2 -46.034225 +/- 0.002539 1.013420 +/- 0.002835 0.0220
qmc.g005 series 2 -46.046731 +/- 0.002963 1.018337 +/- 0.004109 0.0221
qmc.g006 series 2 -46.047133 +/- 0.001958 1.021483 +/- 0.003082 0.0222
qmc.g007 series 2 -45.476146 +/- 0.002065 1.070456 +/- 0.003133 0.0235
The degenerate twists grouped by set are :math:`\{0\}`,
:math:`\{1,2,4\}`, :math:`\{3,5,6\}`, and :math:`\{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:
::
>qmca -a -w '1 3 3 1' -q ev -e 30 *g000*2*sc* *g001*2*sc* *g003*2*sc* *g007*2*sc*
LocalEnergy Variance ratio
avg series 2 -45.873631 +/- 0.001044 1.028769 +/- 0.001520 0.0224
yielding a total energy value of :math:`-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 :ref:`appendix-a`.
.. _qmca-output-units:
Setting output units
~~~~~~~~~~~~~~~~~~~~
Estimates outputted by ``qmca`` are in Hartree units by default. The
output units for energetic quantities can be changed by using the ``-u``
option.
Energy in Hartrees:
::
>qmca -q e -u Ha -e 20 qmc.s002.scalar.dat
qmc series 2 LocalEnergy = -46.032960 +/- 0.002077
Energy in electron volts:
::
>qmca -q e -u eV -e 20 qmc.s002.scalar.dat
qmc series 2 LocalEnergy = -1252.620565 +/- 0.056521
Energy in Rydbergs:
::
>qmca -q e -u rydberg -e 20 qmc.s002.scalar.dat
qmc series 2 LocalEnergy = -92.065919 +/- 0.004154
Energy in kilojoules per mole:
::
>qmca -q e -u kj_mol -e 20 qmc.s002.scalar.dat
qmc series 2 LocalEnergy = -120859.512998 +/- 5.453431
.. _qmca-fast-trace-plot:
Speeding up trace plotting
~~~~~~~~~~~~~~~~~~~~~~~~~~
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
:math:`\mathcal{O}(M^2)`, with :math:`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 -q e -e 20 --noac qmc.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.
.. _qmca-short-example:
Short usage examples
~~~~~~~~~~~~~~~~~~~~
Plotting a trace of the local energy:
::
>qmca -t -q e *scalar*
Applying an equilibration cutoff to VMC data (series 0):
::
>qmca -q e -e 30 *s000.scalar*
Applying the same equilibration cutoff to VMC and DMC data (series 0, 1, 2):
::
>qmca -q e -e 20 *scalar*
Applying different equilibration cutoffs to VMC and DMC data (series 0, 1, 2):
::
>qmca -q e -e '30 20 40' *scalar*
Obtaining the energy, variance, and variance/energy ratio for all series:
::
>qmca -q ev -e 30 *scalar*
Overlaying plots of mean + error bar for energy and variance for separate
two- and three-body Jastrow optimization runs:
::
>qmca -po -q ev ./optJ2/*scalar* ./optJ3/*scalar*
Obtaining the acceptance ratio:
::
>qmca -q ar -e 30 *scalar*
Obtaining the average DMC walker population:
::
>qmca -q nw -e 400 *s002.dmc.dat
Obtaining the MC efficiency:
::
>qmca -q eff -e 30 *scalar*
Obtaining the total wall clock time per series:
::
>qmca -q tt -e 0 *scalar*
Obtaining the average wall clock time spent per block:
::
>qmca -q bc -e 0 *scalar*
Obtaining a subset of desired quantities:
::
>qmca -q 'e v ar eff' -e 30 *scalar*
Obtaining all available quantities:
::
>qmca -e 30 *scalar*
Obtaining the twist-averaged total energy with uniform weights:
::
>qmca -a -q e -e 40 *g*s002.scalar.dat
Obtaining the twist-averaged total energy with specific weights:
::
>qmca -a -w '1 3 3 1' -q e -e 40 *g*s002.scalar.dat
Obtaining the local, kinetic, and potential energies in eV:
::
>qmca -q ekp -e 30 -u eV *scalar*
.. _qmca-production-checklist:
Production quality checklist
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#. 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 -q ev *scalar*``. For
pseudopotential systems without a Jastrow, the variance/energy ratio
should not exceed :math:`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 :math:`0.04` Ha for
pseudopotential systems. A good Jastrow is indicated by a
variance/energy ratio in the range of :math:`0.01-0.03` Ha. A value
less than :math:`0.01` Ha is difficult to achieve.
#. Confirm that the optimization has converged by plotting the energy
and variance vs. optimization series (``qmca -p -q ev *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 :math:`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
:math:`\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 -q nw *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 :math:`0.002-0.01` Ha\ :math:`^{-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
:math:`2 \times 2 \times 2` twist grid in a
:math:`2 \times 2 \times 2` tiled cell is equivalent to a
:math:`1 \times 1 \times 1` twist grid in a
:math:`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).
.. _qmcfit:
Using the qmc-fit tool for statistical time step extrapolation, trial wavefunction optimization and curve fitting
-----------------------------------------------------------------------------------------------------------------
The ``qmc-fit`` tool is used to provide statistical estimates of curve-fitting parameters based on QMCPACK data.
``qmc-fit`` is currently estimates fitting parameters related to time step extrapolation and trial wavefunction
optimization (optimal U for DFT+U, EXX fractions) and supports many types of fitted curves (e.g., Morse
potential binding curves and various equation-of-state fitting curves). An overview of all supported input flags to
``qmc-fit`` can be obtained by typing ``qmc-fit -h`` at the command line:
::
>qmc-fit -h
usage: qmc-fit [-h] [-f FIT_FUNCTION] (-t TIMESTEPS | -u HUBBARDS | --exx EXX | --eos EOS) [-s SERIES_START] [-e EQUILS] [-b REBLOCK_FACTORS] [--noplot] {ts,u,eos} scalar_files [scalar_files ...]
This utility provides a fit to the one-dimensional parameter scans of QMC observables. Currently, the functionality in place is to fit linear/quadratic polynomial fits to the timestep VMC/DMC studies and single parameter optimization of trial wavefunctions using DMC local
energies and quadratic, cubic and quartic fits and equation-of-state and Morse potential fits.
positional arguments:
{ts,u,eos} One dimensional parameter used to fit QMC local energies. Options are ts for timestep and u for hubbard_u parameter fitting
scalar_files Scalar files used in the fit. An explicit list of scalar files with space or a wildcard (e.g. dmc*/dmc.s001.scalar.dat) is acceptable.
optional arguments:
-h, --help show this help message and exit
-f FIT_FUNCTION, --fit FIT_FUNCTION
Fitting function, options are (for each fit type) ts:{linear, quadratic, sqrt} u:{cubic, quadratic, quartic} eos:{birch, morse, murnaghan, vinet}. (default: linear)
-t TIMESTEPS Timesteps corresponding to scalar files, excluding any prior to --series_start (default: None)
-u HUBBARDS Hubbard U values (eV) (default: None)
--exx EXX EXX ratios (default: None)
--eos EOS Structural parameter for EOS fitting (volume or distance) (default: None)
-s SERIES_START, --series_start SERIES_START
Series number for first DMC run. Use to exclude prior VMC scalar files if they have been provided (default: None)
-e EQUILS, --equils EQUILS
Equilibration lengths corresponding to scalar files, excluding any prior to --series_start. Can be a single value for all files. If not provided, equilibration periods will be estimated. (default: None)
-b REBLOCK_FACTORS, --reblock_factors REBLOCK_FACTORS
Reblocking factors corresponding to scalar files, excluding any prior to --series_start. Can be a single value for all files. If not provided, reblocking factors will be estimated. (default: None)
--noplot Do not show plots. (default: False)
The jackknife statistical technique
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 :math:`N` statistical variables :math:`\{x_n\}_{n=1}^N` that
have been outputted by one or more simulation runs. If we have :math:`M`
samples of each of the :math:`N` variables, then the mean values of each
these variables can be estimated in the standard way, that is,
:math:`\bar{x}_n\approx \tfrac{1}{M}\sum_{m=1}^Mx_{nm}`.
Suppose we are interested in :math:`P` statistical quantities
:math:`\{y_p\}_{p=1}^P` that are related to the original :math:`N`
variables by a known multidimensional function :math:`F`:
.. math::
:label: eq46
\begin{aligned}
y_1,y_2,\ldots,y_P &= F(x_1,x_2,\ldots,x_N)\quad \textrm{or} \nonumber \\
\vec{y} &= F(\vec{x})\:.\end{aligned}
The relationship implied by :math:`F` is completely general. For
example, the :math:`\{x_n\}` might be elements of a matrix with
:math:`\{y_p\}` being the eigenvalues, or :math:`F` might be a fitting
procedure for :math:`N` energies at different time steps with :math:`P`
fitting parameters. An approximate guess at the mean value of
:math:`\vec{y}` can be obtained by evaluating :math:`F` at the mean
value of :math:`\vec{x}` (i.e. :math:`F(\bar{x}_1\ldots\bar{x}_N)`), but
with this approach we have no way to estimate the statistical error bar
of any :math:`\bar{y}_p`.
In the jackknife procedure, the statistical variability intrinsic to the
underlying data :math:`\{x_n\}` is used to obtain estimates of the mean
and error bar of :math:`\{y_p\}`. We first construct a new set of
:math:`x` statistical data by taking the average over all samples but
one:
.. math::
:label: eq47
\tilde{x}_{nm} = \frac{1}{N-1}(N\bar{x}_n-x_{nm})\qquad m\in [1,M]\:.
The result is a distribution of approximate :math:`x` mean values. These
are used to construct a distribution of approximate means for :math:`y`:
.. math::
:label: eq48
\tilde{y}_{1m},\ldots,\tilde{y}_{Pm} = F(\tilde{x}_{1m},\ldots,\tilde{x}_{Nm}) \qquad m\in [1,M]\:.
Estimates for the mean and error bar of the quantities of
interest can finally be obtained using the following formulas:
.. math::
:label: eq49
\begin{aligned}
\bar{y}_p &= \frac{1}{M}\sum_{m=1}^M\tilde{y}_{pm}\:. \\
\sigma_{y_p} &= \sqrt{\frac{M-1}{M}\left(\sum_{m=1}^M\tilde{y}_{pm}^2-M\bar{y}_p^2\right)}\:.\end{aligned}
Performing time step extrapolation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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:
:math:`0.04,~0.02,~0.01,~0.005,~0.0025,~0.00125` Ha\ :math:`^{-1}`. For
a typical production pseudopotential study, time steps in the range of
:math:`0.02-0.002` Ha\ :math:`^{-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 -q e -e 0 ./qmc*/*scalar*``), it is clear that
an equilibration period of :math:`30` blocks is reasonable. Approximate
autocorrelation lengths are also obtained with ``qmca``:
::
>qmca -e 30 -q e --sac ./qmc*/qmc.g000.s002.scalar.dat
./qmc_tm_0.00125/qmc.g000 series 2 LocalEnergy = -3848.234513 +/- 0.055754 1.7
./qmc_tm_0.00250/qmc.g000 series 2 LocalEnergy = -3848.237614 +/- 0.055432 2.2
./qmc_tm_0.00500/qmc.g000 series 2 LocalEnergy = -3848.349741 +/- 0.069729 2.8
./qmc_tm_0.01000/qmc.g000 series 2 LocalEnergy = -3848.274596 +/- 0.126407 3.9
./qmc_tm_0.02000/qmc.g000 series 2 LocalEnergy = -3848.539017 +/- 0.075740 2.4
./qmc_tm_0.04000/qmc.g000 series 2 LocalEnergy = -3848.976424 +/- 0.075305 1.8
The autocorrelation must be removed from the data before jackknifing, so
we will reblock the data by a factor of 4.
The ``qmc-fit`` tool can be used in the following way to obtain a linear
time step fit of the data:
::
>qmc-fit ts -e 30 -b 4 -s 2 -t '0.00125 0.0025 0.005 0.01 0.02 0.04' ./qmc*/*scalar*
fit function : linear
fitted formula: (-3848.193 +/- 0.037) + (-18.95 +/- 1.95)*t
intercept : -3848.193 +/- 0.037 Ha
The input arguments are as follows: ``ts`` indicates we are performing a
time step fit, ``-e 30`` is the equilibration period removed from each
set of scalar data, ``-b 4`` 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), ``-s 2`` 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 :math:`-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
:numref:`fig12`.
Different fitting functions are supported via the ``-f`` option.
Currently supported options include ``linear`` (:math:`a+bt`),
``quadratic`` (:math:`a+bt+ct^2`), and ``sqrt``
(:math:`a+b\sqrt{t}+ct`). Results for a quadratic fit are shown
subsequently and in the bottom panel of :numref:`fig12`.
::
>qmc-fit ts -f quadratic -e30 -b4 -s2 -t '0.00125 0.0025 0.005 0.01 0.02 0.04' ./qmc*/*scalar*
fit function : quadratic
fitted formula: (-3848.245 +/- 0.047) + (-7.25 +/- 8.33)*t + (-285.00 +/- 202.39)*t^2
intercept : -3848.245 +/- 0.047 Ha
In this case, we find a zero time step estimate of :math:`-3848.25(5)`
Ha\ :math:`^{-1}`. A time step of :math:`0.04` Ha\ :math:`^{-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 :math:`-3848.22(4)` Ha\ :math:`^{-1}`. Note that quadratic
extrapolation can result in intrinsically larger uncertainty in the
extrapolated value. For example, when the :math:`0.04` Ha\ :math:`^{-1}`
point is excluded, the uncertainty grows by 50% and we obtain an
estimated value of :math:`-3848.28(7)` instead.
.. image:: /figs/qmcfit_timestep_linear.png
:width: 400
:align: center
.. _fig12:
.. figure:: /figs/qmcfit_timestep_quadratic.png
:width: 400
:align: center
Linear (top) and quadratic (bottom) time step fits to DMC data for a 32-atom supercell of MnO obtained with ``qmc-fit``. Zero time step estimates are indicated by the red data point on the left side of either panel.
Performing trial wavefunction optimization fitting, e.g., to find optimal DFT+U
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In this section, we use a 24-atom supercell of monolayer FeCl\ :sub:`2`\ as an example system for wavefunction optimization
fitting. Using single determinant DFT wavefunctions, a practical method to perform wavefunction optimization is done through
scanning the Hubbard-U parameter in a DFT+U calculation used to generate the trial wavefunction. Similarly, one can also scan
different exact exchange ratio parameters in hybrid-DFT calculations. Here, we will show an example of this fitting for the
Hubbard-U parameter, but the same procedure can be applied to any single-parameter scans of trial wavefunctions. Data for this
system has been collected in DMC using the following sequence of Hubbard-U values on Fe-d orbitals: :math:`0 1 2 3 4 5` eV. Some
non-zero U value often minimizes the DMC energy, but optimized U values have limited transferability across different systems.
Similar to the procedure for performing timestep statistical fitting, the quality of the input statistics must be checked using
``qmca`` utility to determine the reblocking factor and equilibration periods. Assuming that an equilibration period of initial 50
steps, ``-e 50``, and a reblocking period of 4, ``-b 6``, is sufficient to remove correlations in the statistical local energies,
the ``qmc-fit`` tool can be used in the following way to obtain a quadratic fit of the data:
::
>qmc-fit u -e 50 -b 6 -u "0 1 2 3 4 5" -f quadratic dmc_u_*/dmc.s001.scalar.dat
fit function : quadratic
fitted formula: (-1230.1071 +/- 0.0045) + (-0.0683 +/- 0.0040)*t + (0.00883 +/- 0.00077)*t^2
root 1 minimum_u : 3.87 +/- 0.14 eV
root 1 minimum_e : -1230.2391 +/- 0.0026 Ha
root 1 curvature : 0.0177 +/- 0.0015
Here, ``qmc-fit u`` indicates we are performing a Hubbard-U/exact-exchange ratio fit,
``-u`` option provides a list of Hubbard-U values :math:`0 1 2 3 4 5` corresponding to the auto-sorted
dmc scalar files with wildcard ``dmc_u_*/dmc.s001.scalar.dat``. Here, ``qmc-fit`` command is invoked at a
directory where folders such as ``dmc_u_0_2x2x1, dmc_u_1_2x2x1, dmc_u_2_2x2x1`` reside.
Here, the text output provides the U value (``minimum_u``) and local energies (``minimum_e``)
at the minima of the polynomial which falls within the range of Hubbard-U values provided in the command
line, e.g. from 0 to 5. Therefore, a U value of :math:`3.8(1)` eV minimizes the DMC energy of the system.
The ``curvature`` is printed for informative purposes only, but a curvature with small error bar could
indicate a higher quality polynomial fit. Similar to the timestep fit, a plot of the fit will also
produced as default where the minima of the polynomial is shown as a red dot as in :numref:`fig13`.
Different fitting functions are supported via the ``-f`` option.
Currently supported options include ``quadratic`` (:math:`a+bt+ct^2`), and
``cubic`` (:math:`a+bt+ct^2+dt^3`) and ``quartic`` (:math:`a+bt+ct^2+dt^3+et^4`).
An example of a cubic fit is given as below:
::
>qmc-fit u -e 50 -b 6 -u "0 1 2 3 4 5" -f cubic dmc_u_*/dmc.s001.scalar.dat
fit function : cubic
fitted formula: (-1230.1087 +/- 0.0045) + (-0.0608 +/- 0.0073)*t + (0.0047 +/- 0.0033)*t^2 + (0.00055 +/- 0.00041)*t^3
root 1 minimum_u : 3.85 +/- 0.11 eV
root 1 minimum_e : -1230.2415 +/- 0.0033 Ha
root 1 curvature : 0.0221 +/- 0.0034
.. _fig13:
.. figure:: /figs/qmcfit_hubbard_quadratic.png
:width: 400
:align: center
Quadratic Hubbard-U fits to DMC data for a 24-atom supercell of monolayer FeCl\ :sub:`2`\ obtained with ``qmc-fit``. DMC local energy minima are indicated by the red data point on the bottom halves of either panel.
Performing equation-of-state and Morse potential binding curve fits
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For a systematic series of statistical data, such as QMC calculations performed at different interatomic distances, or at a series of volumes
for an equation-of-states calculation, it is advised to perform jackknife fitting to determine quantities such as equilibrium distance, volume and
bulk moduli. For interatomic distances and equation of states fits to QMC calculations, ``qmc-fit`` has the capability to perform Morse and Birch, Murnaghan
and Vinet equation-of-state fits. In this example, we determine the equilibrium volume and bulk modulus of C-diamond using a 16 atom supercell using DMC and Murnaghan
equation-of-state fit. For the 16 atom supercell, we uniformly scan over the volumes between :math:`78.16` and :math:`99.62 A^3`. Assuming that all these DMC calculation
folders are located under the same parent folder and ordered from smaller to the large volume (e.g. dmc_78.16, dmc_80.65 ...), the following script can be used to make a
Murnaghan fit to the DMC energies.
::
>qmc-fit eos -e 50 -b 6 --eos "78.16 80.65 83.20 85.80 88.46 91.16 93.92 96.74 99.62" --fit murnaghan dmc_*/dmc.s001.scalar.dat
fit function : murnaghan
fitted formula: E_inf + B/Bp*V*((V_0/V)**Bp/(Bp-1)+1)-V_0*B/(Bp-1)
minimum_x: 89.00 +/- 0.12
e_inf: -91.2659 +/- 0.0012
B: 0.1053 +/- 0.0050
Bp: 0.000189 +/- 0.000011
pressure: -0.00000 +/- 0.00015
Here, the minimum volume is reported as :math:`89.00 \pm 0.12A^3` consistent with the input volume units. Considering that this is a 16-atom cell, the per atom
quantity would be :math:`5.56 \pm 0.01 A^3` per C. Bulk modulus, B, is reported as :math:`0.1053 \pm 0.005 Ha/A^3`. In SI units, this bulk modulus value corresponds to
:math:`459 \pm 21` GPa. Different fitting functions are supported via the ``-f`` option. Currently supported options include ``Vinet`` , ``Murnaghan``, ``Birch`` and ``Morse``.
For more information and default options, please refer to ``qmc-fit -h``.
.. _fig14:
.. figure:: /figs/qmcfit_eos.png
:width: 400
:align: center
Murnaghan equation-of-state fits to DMC data for a 16-atom supercell of C-diamond obtained with ``qmc-fit``. DMC structural minimum is indicated by the red data point with an error bar smaller than its marker size.
.. _qdens:
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:
.. code-block::
>qdens
Usage: qdens [options] [file(s)]
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).
-f FORMATS, --formats=FORMATS
Format or list of formats for density file output.
Options: dat, xsf, chgcar (default=None).
-e EQUILIBRATION, --equilibration=EQUILIBRATION
Equilibration length in blocks (default=0).
-r REBLOCK, --reblock=REBLOCK
Block coarsening factor; use estimated autocorrelation
length (default=None).
-a, --average Average over files in each series (default=False).
-w WEIGHTS, --weights=WEIGHTS
List of weights for averaging (default=None).
-i INPUT, --input=INPUT
QMCPACK input file containing structure and grid
information (default=None).
-s STRUCTURE, --structure=STRUCTURE
File containing atomic structure (default=None).
-g GRID, --grid=GRID Density grid dimensions (default=None).
-c CELL, --cell=CELL Simulation cell axes (default=None).
--density_cell=DENSITY_CELL
Density cell axes (default=None).
--density_corner=DENSITY_CORNER
Density cell corner (default=None).
--lineplot=LINEPLOT Produce a line plot along the selected dimension: 0,
1, or 2 (default=None).
--noplot Do not show plots interactively (default=False).
--twist_info=TWIST_INFO
Use twist weights in twist_info.dat files or not.
Options: "use", "ignore", "require". "use" means use
when present, "ignore" means do not use, "require"
means must be used (default=use).
Usage examples
~~~~~~~~~~~~~~
Process a single file, excluding the first 40 blocks, and produce XSF files:
::
qdens -v -e 40 -f xsf -i qmc.in.xml qmc.s000.stat.h5
Process files for all available series:
::
qdens -v -e 40 -f xsf -i qmc.in.xml *stat.h5
Combine groups of 10 adjacent statistical blocks together (appropriate if the
estimated autocorrelation time is about 10 blocks):
::
qdens -v -e 40 -r 10 -f xsf -i qmc.in.xml qmc.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``):
::
qdens -v -e '20 20 40' -r '4 4 8' -f xsf -i qmc.in.xml *stat.h5
Produce twist averaged densities (also works with multiple series and reblocking):
::
qdens -v -a -e 40 -f xsf -i qmc.g000.twistnum_0.in.xml qmc.g*.s000.stat.h5
Twist averaging with arbitrary weights can be performed via the ``-w``
option in a fashion identical to ``qmca``.
Files produced
~~~~~~~~~~~~~~
Look for files with names and extensions similar to:
::
qmc.s000.SpinDensity_u.xsf
qmc.s000.SpinDensity_u-err.xsf
qmc.s000.SpinDensity_u+err.xsf
qmc.s000.SpinDensity_d.xsf
qmc.s000.SpinDensity_d-err.xsf
qmc.s000.SpinDensity_d+err.xsf
qmc.s000.SpinDensity_u+d.xsf
qmc.s000.SpinDensity_u+d-err.xsf
qmc.s000.SpinDensity_u+d+err.xsf
qmc.s000.SpinDensity_u-d.xsf
qmc.s000.SpinDensity_u-d-err.xsf
qmc.s000.SpinDensity_u-d+err.xsf
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.
.. _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
Reference in New Issue View Git Blame Copy Permalink