Hamiltonian and Observables¶
QMCPACK is capable of the simultaneous measurement of the Hamiltonian and many other quantum operators. The Hamiltonian attains a special status among the available operators (also referred to as observables) because it ultimately generates all available information regarding the quantum system. This is evident from an algorithmic standpoint as well since the Hamiltonian (embodied in the projector) generates the imaginary time dynamics of the walkers in DMC and reptation Monte Carlo (RMC).
This section covers how the Hamiltonian can be specified, component by component, by the user in the XML format native to qmcpack. It also covers the input structure of statistical estimators corresponding to quantum observables such as the density, static structure factor, and forces.
The Hamiltonian¶
The manybody Hamiltonian in Hartree units is given by
Here, the sums indexed by \(i/j\) are over quantum particles, while \(\ell/m\) are reserved for classical particles. Often the quantum particles are electrons, and the classical particles are ions, though is not limited in this way. The mass of each quantum particle is denoted \(m_i\), \(v^{qq}/v^{qc}/v^{cc}\) are pair potentials between quantumquantum/quantumclassical/classicalclassical particles, and \(v^{ext}\) denotes a purely external potential.
QMCPACK is designed modularly so that any potential can be supported with minimal additions to the code base. Potentials currently supported include Coulomb interactions in open and periodic boundary conditions, the MPC potential, nonlocal pseudopotentials, helium pair potentials, and various model potentials such as hard sphere, Gaussian, and modified PoschlTeller.
Reference information and examples for the <hamiltonian/>
XML
element are provided subsequently. Detailed descriptions of the input
for individual potentials is given in the sections that follow.
hamiltonian
element:
parent elements:
simulation, qmcsystem
child elements:
pairpot extpot estimator constant
(deprecated)
attributes:
Name
Datatype
Values
Default
Description
name/id
\(^o\)text
anything
h0
Unique id for this Hamiltonian instance
type
\(^o\)text
generic
No current function
role
\(^o\)text
primary/extra
extra
Designate as Hamiltonian or not
source
\(^o\)text
particleset.name
i
Identify classical
particleset
target
\(^o\)text
particleset.name
e
Identify quantum
particleset
default
\(^o\)boolean
yes/no
yes
Include kinetic energy term implicitly
Additional information:
target: Must be set to the name of the quantum
particleset
. The default value is typically sufficient. In normal usage, no other attributes are provided.
<hamiltonian target="e">
<pairpot name="ElecElec" type="coulomb" source="e" target="e"/>
<pairpot name="ElecIon" type="coulomb" source="i" target="e"/>
<pairpot name="IonIon" type="coulomb" source="i" target="i"/>
</hamiltonian>
<hamiltonian target="e">
<pairpot name="ElecElec" type="coulomb" source="e" target="e"/>
<pairpot name="PseudoPot" type="pseudo" source="i" wavefunction="psi0" format="xml">
<pseudo elementType="Li" href="Li.xml"/>
<pseudo elementType="H" href="H.xml"/>
</pairpot>
<pairpot name="IonIon" type="coulomb" source="i" target="i"/>
</hamiltonian>
Pair potentials¶
Many pair potentials are supported. Though only the most commonly used pair potentials are covered in detail in this section, all currently available potentials are listed subsequently. If a potential you desire is not listed, or is not present at all, feel free to contact the developers.
pairpot
factory element:
parent elements:
hamiltonian
child elements:
type
attribute
type options
coulomb
Coulomb/Ewald potential
pseudo
Semilocal pseudopotential
mpc
Model periodic Coulomb interaction/correction
cpp
Core polarization potential
skpot
Unknown
shared attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
See above
0
Select pairpot type
name
\(^r\)text
Anything
any
Unique name for this pairpot
source
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
target
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
units
\(^o\)text
hartree
No current function
Additional information:
type: Used to select the desired pair potential. Must be selected from the list of type options.
name: A unique name used to identify this pair potential. Block averaged output data will appear under this name in
scalar.dat
and/orstat.h5
files.source/target: These specify the particles involved in a pair interaction. If an interaction is between classical (e.g., ions) and quantum (e.g., electrons),
source
/target
should be the name of the classical/quantumparticleset
.Only
Coulomb, pseudo
, andmpc
are described in detail in the following subsections. The older or lessused types (cpp, skpot
) are not covered.Available only if
QMC_CUDA
is not defined:skpot
.Available only if
OHMMS_DIM==3
:mpc, vhxc, pseudo
.Available only if
OHMMS_DIM==3
andQMC_CUDA
is not defined:cpp
.
Coulomb potentials¶
The bare Coulomb potential is used in open boundary conditions:
When periodic boundary conditions are selected, Ewald summation is used automatically:
The sum indexed by \(L\) is over all nonzero simulation cell lattice vectors. In practice, the Ewald sum is broken into short and longrange parts in a manner optimized for efficiency (see [NC95]) for details.
For information on how to set the boundary conditions, consult Specifying the system to be simulated.
pairpot type=coulomb
element:
parent elements:
hamiltonian
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
coulomb
Must be coulomb
name/id
\(^r\)text
anything
ElecElec
Unique name for interaction
source
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
target
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
pbc
\(^o\)boolean
yes/no
yes
Use Ewald summation
physical
\(^o\)boolean
yes/no
yes
Hamiltonian(yes)/Observable(no)
gpu
boolean
yes/no
depend
Offload computation to GPU
forces
boolean
yes/no
no
Deprecated
Additional information:
type/source/target: See description for the previous generic
pairpot
factory element.name: Traditional userspecified names for electronelectron, electronion, and ionion terms are
ElecElec
,ElecIon
, andIonIon
, respectively. Although any choice can be used, the data analysis tools expect to find columns in*.scalar.dat
with these names.pbc: Ewald summation will not be performed if
simulationcell.bconds== n n n
, regardless of the value ofpbc
. Similarly, thepbc
attribute can only be used to turn off Ewald summation ifsimulationcell.bconds!= n n n
. The default value is recommended.physical: If
physical==yes
, this pair potential is included in the Hamiltonian and will factor into theLocalEnergy
reported by QMCPACK and also in the DMC branching weight. Ifphysical==no
, then the pair potential is treated as a passive observable but not as part of the Hamiltonian itself. As such it does not contribute to the outputtedLocalEnergy
. Regardless of the value ofphysical
output data will appear inscalar.dat
in a column headed byname
.gpu: When not specified, use the
gpu
attribute ofparticleset
.
<pairpot name="ElecElec" type="coulomb" source="e" target="e"/>
<pairpot name="ElecIon" type="coulomb" source="i" target="e"/>
<pairpot name="IonIon" type="coulomb" source="i" target="i"/>
Pseudopotentials¶
QMCPACK supports pseudopotentials in semilocal form, which is local in the radial coordinate and nonlocal in angular coordinates. When all angular momentum channels above a certain threshold (\(\ell_{max}\)) are well approximated by the same potential (\(V_{\bar{\ell}}\equiv V_{loc}\)), the pseudopotential separates into a fully local channel and an angularly nonlocal component:
Here the electron/ion index is \(i/j\), and only one type of ion is shown for simplicity.
Evaluation of the localized pseudopotential energy \(\Psi_T^{1}V^{PP}\Psi_T\) requires additional angular integrals. These integrals are evaluated on a randomly shifted angular grid. The size of this grid is determined by \(\ell_{max}\). See [MSC91] for further detail.
uses the FSAtom pseudopotential file format associated with the “Free
Software Project for Atomicscale Simulations” initiated in 2002. See
http://www.tddft.org/fsatom/manifest.php for more information. The
FSAtom format uses XML for structured data. Files in this format do not
use a specific identifying file extension; instead they are simply
suffixed with “.xml
.” The tabular data format of CASINO is also
supported.
In addition to the semilocal pseudopotential above, spinorbit interactions can also be included through the use of spinorbit pseudopotentials. The spinorbit contribution can be written as
Here, \(\vec{s}\) is the spin operator. For each atom with a spinorbit contribution,
the radial functions \(V_{\ell}^{\rm SO}\) can be included in the pseudopotential
“.xml
” file.
pairpot type=pseudo
element:
parent elements:
hamiltonian
child elements:
pseudo
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
pseudo
Must be pseudo
name/id
\(^r\)text
anything
PseudoPot
No current function
source
\(^r\)text
particleset.name
i
Ion
particleset
name
target
\(^r\)text
particleset.name
hamiltonian.target
Electron
particleset
name
pbc
\(^o\)boolean
yes/no
yes*
Use Ewald summation
forces
boolean
yes/no
no
Deprecated
wavefunction
\(^r\)text
wavefunction.name
invalid
Identify wavefunction
format
\(^r\)text
xml/table
table
Select file format
algorithm
\(^o\)text
batched/nonbatched
batched
Choose NLPP algorithm
DLA
\(^o\)text
yes/no
no
Use determinant localization approximation
physicalSO
\(^o\)boolean
yes/no
yes
Include the SO contribution in the local energy
Additional information:
type/source/target See description for the generic
pairpot
factory element.name: Ignored. Instead, default names will be present in
*scalar.dat
output files when pseudopotentials are used. The fieldLocalECP
refers to the local part of the pseudopotential. If nonlocal channels are present, aNonLocalECP
field will be added that contains the nonlocal energy summed over all angular momentum channels.pbc: Ewald summation will not be performed if
simulationcell.bconds== n n n
, regardless of the value ofpbc
. Similarly, thepbc
attribute can only be used to turn off Ewald summation ifsimulationcell.bconds!= n n n
.format: If
format
==table, QMCPACK looks for*.psf
files containing pseudopotential data in a tabular format. The files must be named after the ionic species provided inparticleset
(e.g.,Li.psf
andH.psf
). Ifformat
==xml, additionalpseudo
child XML elements must be provided (see the following). These elements specify individual file names and formats (both the FSAtom XML and CASINO tabular data formats are supported).algorithm The
nonbatched
algorithm evaluates the ratios of wavefunction components together for each quadrature point and then one point after another. Thebatched
algorithm evaluates the ratios of quadrature points together for each wavefunction component and then one component after another. Internally, it usesVirtualParticleSet
for quadrature points. Hybrid orbital representation has an extra optimization enabled when using the batched algorithm. When OpenMP offload build is enabled, the default value isbatched
. Otherwise,nonbatched
is the default.DLA Determinant localization approximation (DLA) [ZBMAlfe19] uses only the fermionic part of the wavefunction when calculating NLPP.
physicalSO If the spinorbit components are included in the
.xml
file, this flag allows control over whether the SO contribution is included in the local energy.
<pairpot name="PseudoPot" type="pseudo" source="i" wavefunction="psi0" format="psf"/>
<pairpot name="PseudoPot" type="pseudo" source="i" wavefunction="psi0" format="xml">
<pseudo elementType="Li" href="Li.xml"/>
<pseudo elementType="H" href="H.xml"/>
</pairpot>
<pairpot name="PseudoPot" type="pseudo" source="i" wavefunction="psi0" format="xml" physicalSO="no">
<pseudo elementType="Pb" href="Pb.xml"/>
</pairpot>
Details of <pseudo/>
input elements are shown in the following. It
is possible to include (or construct) a full pseudopotential directly in
the input file without providing an external file via href
. The full
XML format for pseudopotentials is not yet covered.
pseudo
element:
parent elements:
pairpot type=pseudo
child elements:
header local grid
attributes:
Name
Datatype
Values
Default
Description
elementType/symbol
\(^r\)text
groupe.name
none
Identify ionic species
href
\(^r\)text
filepath
none
Pseudopotential file path
format
\(^r\)text
xml/casino
xml
Specify file format
cutoff
\(^o\)real
Nonlocal cutoff radius
lmax
\(^o\)integer
Largest angular momentum
nrule
\(^o\)integer
Integration grid order
llocal
\(^o\)integer
Override local channel
<pseudo elementType="Li" href="Li.xml"/>
MPC Interaction/correction¶
The MPC interaction is an alternative to direct Ewald summation. The MPC
corrects the exchange correlation hole to more closely match its
thermodynamic limit. Because of this, the MPC exhibits smaller
finitesize errors than the bare Ewald interaction, though a few
alternative and competitive finitesize correction schemes now exist.
The MPC is itself often used just as a finitesize correction in
postprocessing (set physical=false
in the input).
pairpot type=mpc
element:
parent elements:
hamiltonian
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
mpc
Must be MPC
name/id
\(^r\)text
anything
MPC
Unique name for interaction
source
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
target
\(^r\)text
particleset.name
hamiltonian.target
Identify interacting particles
physical
\(^o\)boolean
yes/no
no
Hamiltonian(yes)/observable(no)
cutoff
real
\(>0\)
30.0
Kinetic energy cutoff
Remarks:
physical
: Typically set tono
, meaning the standard Ewald interaction will be used during sampling and MPC will be measured as an observable for finitesize postcorrection. Ifphysical
isyes
, the MPC interaction will be used during sampling. In this case an electronelectron Coulombpairpot
element should not be supplied.Developer note: Currently the
name
attribute for the MPC interaction is ignored. The name is always reset toMPC
.
<pairpot type="MPC" name="MPC" source="e" target="e" ecut="60.0" physical="no"/>
General estimators¶
A broad range of estimators for physical observables are available in QMCPACK.
The following sections contain input details for the total number
density (density
), number density resolved by particle spin
(spindensity
), spherically averaged pair correlation function
(gofr
), static structure factor (sk
), static structure factor
(skall
), energy density (energydensity
), one body reduced
density matrix (dm1b
), \(S(k)\) based kinetic energy correction
(chiesa
), forward walking (ForwardWalking
), and force
(Force
) estimators. Other estimators are not yet covered.
When an <estimator/>
element appears in <hamiltonian/>
, it is
evaluated for all applicable chained QMC runs (e.g.,
VMC\(\rightarrow\)DMC\(\rightarrow\)DMC). Estimators are
generally not accumulated during wavefunction optimization sections. If
an <estimator/>
element is instead provided in a particular
<qmc/>
element, that estimator is only evaluated for that specific
section (e.g., during VMC only).
estimator
factory element:
parent elements:
hamiltonian, qmc
type selector:
type
attribute
type options
density
Density on a grid
spindensity
Spin density on a grid
gofr
Pair correlation function (quantum species)
sk
Static structure factor
SkAll
Static structure factor needed for finite size correction
structurefactor
Species resolved structure factor
species kinetic
Species resolved kinetic energy
latticedeviation
Spatial deviation between two particlesets
momentum
Momentum distribution
energydensity
Energy density on uniform or Voronoi grid
dm1b
One body density matrix in arbitrary basis
chiesa
ChiesaCeperleyMartinHolzmann kinetic energy correction
Force
Family of “force” estimators (see ChiesaCeperleyZhang Force Estimators)
ForwardWalking
Forward walking values for existing estimators
orbitalimages
Create image files for orbitals, then exit
flux
Checks sampling of kinetic energy
localmoment
Atomic spin polarization within cutoff radius
Pressure
No current function
shared attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
See above
0
Select estimator type
name
\(^r\)text
anything
any
Unique name for this estimator
ChiesaCeperleyMartinHolzmann kinetic energy correction¶
This estimator calculates a finitesize correction to the kinetic energy following the formalism laid out in [CCMH06]. The total energy can be corrected for finitesize effects by using this estimator in conjunction with the MPC correction.
estimator type=chiesa
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
chiesa
Must be chiesa
name
\(^o\)text
anything
KEcorr
Always reset to KEcorr
source
\(^o\)text
particleset.name
e
Identify quantum particles
psi
\(^o\)text
wavefunction.name
psi0
Identify wavefunction
<estimator name="KEcorr" type="chiesa" source="e" psi="psi0"/>
Density estimator¶
The particle number density operator is given by
The density
estimator accumulates the number density on a uniform
histogram grid over the simulation cell. The value obtained for a grid
cell \(c\) with volume \(\Omega_c\) is then the average number
of particles in that cell:
estimator type=density
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
density
Must be density
name
\(^r\)text
anything
any
Unique name for estimator
delta
\(^o\)real array(3)
\(0\le v_i \le 1\)
0.1 0.1 0.1
Grid cell spacing, unit coords
x_min
\(^o\)real
\(>0\)
0
Grid starting point in x (Bohr)
x_max
\(^o\)real
\(>0\)
\(\)
lattice[0]
\(\)Grid ending point in x (Bohr)
y_min
\(^o\)real
\(>0\)
0
Grid starting point in y (Bohr)
y_max
\(^o\)real
\(>0\)
\(\)
lattice[1]
\(\)Grid ending point in y (Bohr)
z_min
\(^o\)real
\(>0\)
0
Grid starting point in z (Bohr)
z_max
\(^o\)real
\(>0\)
\(\)
lattice[2]
\(\)Grid ending point in z (Bohr)
potential
\(^o\)boolean
yes/no
no
Accumulate local potential, deprecated
debug
\(^o\)boolean
yes/no
no
No current function
Additional information:
name
: The name provided will be used as a label in thestat.h5
file for the blocked output data. Postprocessing tools expectname="Density."
delta
: This sets the histogram grid size used to accumulate the density:delta="0.1 0.1 0.05"
\(\rightarrow 10\times 10\times 20\) grid,delta="0.01 0.01 0.01"
\(\rightarrow 100\times 100\times 100\) grid. The density grid is written to astat.h5
file at the end of each MC block. If you request many \(blocks\) in a<qmc/>
element, or select a large grid, the resultingstat.h5
file could be many gigabytes in size.*_min/*_max
: Can be used to select a subset of the simulation cell for the density histogram grid. For example if a (cubic) simulation cell is 20 Bohr on a side, setting*_min=5.0
and*_max=15.0
will result in a density histogram grid spanning a \(10\times 10\times 10\) Bohr cube about the center of the box. Use ofx_min, x_max, y_min, y_max, z_min, z_max
is only appropriate for orthorhombic simulation cells with open boundary conditions.When open boundary conditions are used, a
<simulationcell/>
element must be explicitly provided as the first subelement of<qmcsystem/>
for the density estimator to work. In this case the molecule should be centered around the middle of the simulation cell (\(L/2\)) and not the origin (\(0\) since the space within the cell, and hence the density grid, is defined from \(0\) to \(L\)).
<estimator name="Density" type="density" delta="0.05 0.05 0.05"/>
Spin density estimator¶
The spin density is similar to the total density described previously. In this case, the sum over particles is performed independently for each spin component.
estimator type=spindensity
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
spindensity
Must be spindensity
name
\(^r\)text
anything
any
Unique name for estimator
report
\(^o\)boolean
yes/no
no
Write setup details to stdout
parameters:
Name
Datatype
Values
Default
Description
grid
\(^o\)integer array(3)
\(v_i>\)
Grid cell count
dr
\(^o\)real array(3)
\(v_i>\)
Grid cell spacing (Bohr)
cell
\(^o\)real array(3,3)
anything
Volume grid exists in
corner
\(^o\)real array(3)
anything
Volume corner location
center
\(^o\)real array (3)
anything
Volume center/origin location
voronoi
\(^o\)text
particleset.name
Under development
test_moves
\(^o\)integer
\(>=0\)
0
Test estimator with random moves
Additional information:
name
: The name provided will be used as a label in thestat.h5
file for the blocked output data. Postprocessing tools expectname="SpinDensity."
grid
: The grid sets the dimension of the histogram grid. Input like<parameter name="grid"> 40 40 40 </parameter>
requests a \(40 \times 40\times 40\) grid. The shape of individual grid cells is commensurate with the supercell shape.dr
: Thedr
sets the realspace dimensions of grid cell edges (Bohr units). Input like<parameter name="dr"> 0.5 0.5 0.5 </parameter>
in a supercell with axes of length 10 Bohr each (but of arbitrary shape) will produce a \(20\times 20\times 20\) grid. The inputteddr
values are rounded to produce an integer number of grid cells along each supercell axis. Eithergrid
ordr
must be provided, but not both.cell
: Whencell
is provided, a userdefined grid volume is used instead of the global supercell. This must be provided if open boundary conditions are used. Additionally, ifcell
is provided, the user must specify where the volume is located in space in addition to its size/shape (cell
) using either thecorner
orcenter
parameters.corner
: The grid volume is defined as \(corner+\sum_{d=1}^3u_dcell_d\) with \(0<u_d<1\) (“cell” refers to either the supercell or userprovided cell).center
: The grid volume is defined as \(center+\sum_{d=1}^3u_dcell_d\) with \(1/2<u_d<1/2\) (“cell” refers to either the supercell or userprovided cell).corner/center
can be used to shift the grid even ifcell
is not specified. Simultaneous use ofcorner
andcenter
will cause QMCPACK to abort.
<estimator type="spindensity" name="SpinDensity" report="yes">
<parameter name="grid"> 40 40 40 </parameter>
</estimator>
<estimator type="spindensity" name="SpinDensity" report="yes">
<parameter name="grid">
20 20 20
</parameter>
<parameter name="center">
0.0 0.0 0.0
</parameter>
<parameter name="cell">
10.0 0.0 0.0
0.0 10.0 0.0
0.0 0.0 10.0
</parameter>
</estimator>
Pair correlation function, \(g(r)\)¶
The functional form of the speciesresolved radial pair correlation function operator is
where \(N_s\) is the number of particles of species \(s\) and \(V\) is the supercell volume. If \(s=s'\), then the sum is restricted so that \(i_s\ne j_s\).
In QMCPACK, an estimate of \(g_{ss'}(r)\) is obtained as a radial histogram with a set of \(N_b\) uniform bins of width \(\delta r\). This can be expressed analytically as
where the radial coordinate \(r\) is restricted to reside at the bin centers, \(\delta r/2, 3 \delta r/2, 5 \delta r/2, \ldots\).
estimator type=gofr
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
gofr
Must be gofr
name
\(^o\)text
anything
any
No current function
num_bin
\(^r\)integer
\(>1\)
20
# of histogram bins
rmax
\(^o\)real
\(>0\)
10
Histogram extent (Bohr)
dr
\(^o\)real
\(0\)
0.5
No current function
debug
\(^o\)boolean
yes/no
no
No current function
target
\(^o\)text
particleset.name
hamiltonian.target
Quantum particles
source/sources
\(^o\)text array
particleset.name
hamiltonian.target
Classical particles
Additional information:
num_bin:
This is the number of bins in each species pair radial histogram.rmax:
This is the maximum pair distance included in the histogram. The uniform bin width is \(\delta r=\texttt{rmax/num\_bin}\). If periodic boundary conditions are used for any dimension of the simulation cell, then the default value ofrmax
is the simulation cell radius instead of 10 Bohr. For open boundary conditions, the volume (\(V\)) used is 1.0 Bohr\(^3\).source/sources:
If unspecified, only pair correlations between each species of quantum particle will be measured. For each classical particleset specified bysource/sources
, additional pair correlations between each quantum and classical species will be measured. Typically there is only one classical particleset (e.g.,source="ion0"
), but there can be several in principle (e.g.,sources="ion0 ion1 ion2"
).target:
The default value is the preferred usage (i.e.,target
does not need to be provided).Data is output to the
stat.h5
for each QMC subrun. Individual histograms are named according to the quantum particleset and index of the pair. For example, if the quantum particleset is named “e” and there are two species (up and down electrons, say), then there will be three sets of histogram data in eachstat.h5
file namedgofr_e_0_0
,gofr_e_0_1
, andgofr_e_1_1
for upup, updown, and downdown correlations, respectively.
<estimator type="gofr" name="gofr" num_bin="200" rmax="3.0" />
<estimator type="gofr" name="gofr" num_bin="200" rmax="3.0" source="ion0" />
Static structure factor, \(S(k)\)¶
Let
\(\rho^e_{\mathbf{k}}=\sum_j e^{i \mathbf{k}\cdot\mathbf{r}_j^e}\)
be the Fourier space electron density, with \(\mathbf{r}^e_j\) being
the coordinate of the jth electron. \(\mathbf{k}\) is a wavevector
commensurate with the simulation cell. QMCPACK allows the user to
accumulate the static electron structure factor \(S(\mathbf{k})\) at
all commensurate \(\mathbf{k}\) such that
\(\mathbf{k} \leq (LR\_DIM\_CUTOFF) r_c\). \(N^e\) is the
number of electrons, LR_DIM_CUTOFF
is the optimized breakup
parameter, and \(r_c\) is the WignerSeitz radius. It is defined as
follows:
estimator type=sk
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
sk
Must sk
name
\(^r\)text
anything
any
Unique name for estimator
hdf5
\(^o\)boolean
yes/no
no
Output to
stat.h5
(yes) orscalar.dat
(no)
Additional information:
name:
This is the unique name for estimator instance. A data structure of the same name will appear instat.h5
output files.hdf5:
Ifhdf5==yes
, output data for \(S(k)\) is directed to thestat.h5
file (recommended usage). Ifhdf5==no
, the data is instead routed to thescalar.dat
file, resulting in many columns of data with headings prefixed byname
and postfixed by the kpoint index (e.g.,sk_0 sk_1 …sk_1037 …
).This estimator only works in periodic boundary conditions. Its presence in the input file is ignored otherwise.
This is not a speciesresolved structure factor. Additionally, for \(\mathbf{k}\) vectors commensurate with the unit cell, \(S(\mathbf{k})\) will include contributions from the static electronic density, thus meaning it will not accurately measure the electronelectron density response.
<estimator type="sk" name="sk" hdf5="yes"/>
Static structure factor, SkAll
¶
In order to compute the finite size correction to the potential energy,
records of \(\rho(\mathbf{k})\) is required. What sets SkAll
apart from sk
is that SkAll
records \(\rho(\mathbf{k})\) in
addition to \(s(\mathbf{k})\).
estimator type=SkAll
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
sk
Must be sk
name
\(^r\)text
anything
any
Unique name for estimator
source
\(^r\)text
Ion ParticleSet name
None

target
\(^r\)text
Electron ParticleSet name
None

hdf5
\(^o\)boolean
yes/no
no
Output to
stat.h5
(yes) orscalar.dat
(no)
writeionion
\(^o\)boolean
yes/no
no
Writes file rhok_IonIon.dat containing \(s(\mathbf{k})\) for the ions
Additional information:
name:
This is the unique name for estimator instance. A data structure of the same name will appear instat.h5
output files.hdf5:
Ifhdf5==yes
, output data is directed to thestat.h5
file (recommended usage). Ifhdf5==no
, the data is instead routed to thescalar.dat
file, resulting in many columns of data with headings prefixed byrhok
and postfixed by the kpoint index.This estimator only works in periodic boundary conditions. Its presence in the input file is ignored otherwise.
This is not a speciesresolved structure factor. Additionally, for \(\mathbf{k}\) vectors commensurate with the unit cell, \(S(\mathbf{k})\) will include contributions from the static electronic density, thus meaning it wil not accurately measure the electronelectron density response.
<estimator type="skall" name="SkAll" source="ion0" target="e" hdf5="yes"/>
Species kinetic energy¶
Record speciesresolved kinetic energy instead of the total kinetic
energy in the Kinetic
column of scalar.dat. SpeciesKineticEnergy
is arguably the simplest estimator in QMCPACK. The implementation of
this estimator is detailed in
manual/estimator/estimator_implementation.pdf
.
estimator type=specieskinetic
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
specieskinetic
Must be specieskinetic
name
\(^r\)text
anything
any
Unique name for estimator
hdf5
\(^o\)boolean
yes/no
no
Output to
stat.h5
(yes)
<estimator type="specieskinetic" name="skinetic" hdf5="no"/>
Lattice deviation estimator¶
Record deviation of a group of particles in one particle set (target) from a group of particles in another particle set (source).
estimator type=latticedeviation
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
latticedeviation
Must be latticedeviation
name
\(^r\)text
anything
any
Unique name for estimator
hdf5
\(^o\)boolean
yes/no
no
Output to
stat.h5
(yes)
per_xyz
\(^o\)boolean
yes/no
no
Directionally resolved (yes)
source
\(^r\)text
e/ion0/…
no
source particleset
sgroup
\(^r\)text
u/d/…
no
source particle group
target
\(^r\)text
e/ion0/…
no
target particleset
tgroup
\(^r\)text
u/d/…
no
target particle group
Additional information:
source
: The “reference” particleset to measure distances from; actual reference points are determined together withsgroup
.sgroup
: The “reference” particle group to measure distances from.source
: The “target” particleset to measure distances to.sgroup
: The “target” particle group to measure distances to. For example, in Listing 32 the distance from the up electron (“u”) to the origin of the coordinate system is recorded.per_xyz
: Used to record directionresolved distance. In Listing 32, the x,y,z coordinates of the up electron will be recorded separately ifper_xyz=yes
.hdf5
: Used to record particleresolved distances in the h5 file ifgdf5=yes
.
<particleset name="e" random="yes">
<group name="u" size="1" mass="1.0">
<parameter name="charge" > 1 </parameter>
<parameter name="mass" > 1.0 </parameter>
</group>
<group name="d" size="1" mass="1.0">
<parameter name="charge" > 1 </parameter>
<parameter name="mass" > 1.0 </parameter>
</group>
</particleset>
<particleset name="wf_center">
<group name="origin" size="1">
<attrib name="position" datatype="posArray" condition="0">
0.00000000 0.00000000 0.00000000
</attrib>
</group>
</particleset>
<estimator type="latticedeviation" name="latdev" hdf5="yes" per_xyz="yes"
source="wf_center" sgroup="origin" target="e" tgroup="u"/>
Energy density estimator¶
An energy density operator, \(\hat{\mathcal{E}}_r\), satisfies
where the integral is over all space and \(\hat{H}\) is the Hamiltonian. In QMCPACK, the energy density is split into kinetic and potential components
with each component given by
Here, \(r_i\) and \(\tilde{r}_\ell\) represent electron and ion positions, respectively; \(\hat{p}_i\) is a single electron momentum operator; and \(\hat{v}^{ee}(r_i,r_j)\), \(\hat{v}^{eI}(r_i,\tilde{r}_\ell)\), and \(\hat{v}^{II}(\tilde{r}_\ell,\tilde{r}_m)\) are the electronelectron, electronion, and ionion pair potential operators (including nonlocal pseudopotentials, if present). This form of the energy density is size consistent; that is, the partially integrated energy density operators of wellseparated atoms gives the isolated Hamiltonians of the respective atoms. For periodic systems with twistaveraged boundary conditions, the energy density is formally correct only for either a set of supercell kpoints that correspond to realvalued wavefunctions or a kpoint set that has inversion symmetry around a kpoint having a realvalued wavefunction. For more information about the energy density, see [KYKC13].
In QMCPACK, the energy density can be accumulated on piecewise uniform 3D grids in generalized Cartesian, cylindrical, or spherical coordinates. The energy density integrated within Voronoi volumes centered on ion positions is also available. The total particle number density is also accumulated on the same grids by the energy density estimator for convenience so that related quantities, such as the regional energy per particle, can be computed easily.
estimator type=EnergyDensity
element:
parent elements:
hamiltonian, qmc
child elements:
reference_points, spacegrid
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
EnergyDensity
Must be EnergyDensity
name
\(^r\)text
anything
Unique name for estimator
dynamic
\(^r\)text
particleset.name
Identify electrons
static
\(^o\)text
particleset.name
Identify ions
ion_points
\(^o\)text
yes/no
no
Separate ion energy density onto point field
Additional information:
name:
Must be unique. A dataset with blocked statistical data for the energy density will appear in thestat.h5
files labeled asname
.Important: in order for the estimator to work, a traces XML input element (<traces array=”yes” write=”no”/>) must appear following the <qmcsystem/> element and prior to any <qmc/> element.
<estimator type="EnergyDensity" name="EDcell" dynamic="e" static="ion0">
<spacegrid coord="cartesian">
<origin p1="zero"/>
<axis p1="a1" scale=".5" label="x" grid="1 (.05) 1"/>
<axis p1="a2" scale=".5" label="y" grid="1 (.1) 1"/>
<axis p1="a3" scale=".5" label="z" grid="1 (.1) 1"/>
</spacegrid>
</estimator>
<estimator type="EnergyDensity" name="EDatom" dynamic="e" static="ion0">
<reference_points coord="cartesian">
r1 1 0 0
r2 0 1 0
r3 0 0 1
</reference_points>
<spacegrid coord="spherical">
<origin p1="ion01"/>
<axis p1="r1" scale="6.9" label="r" grid="0 1"/>
<axis p1="r2" scale="6.9" label="phi" grid="0 1"/>
<axis p1="r3" scale="6.9" label="theta" grid="0 1"/>
</spacegrid>
<spacegrid coord="spherical">
<origin p1="ion02"/>
<axis p1="r1" scale="6.9" label="r" grid="0 1"/>
<axis p1="r2" scale="6.9" label="phi" grid="0 1"/>
<axis p1="r3" scale="6.9" label="theta" grid="0 1"/>
</spacegrid>
</estimator>
<estimator type="EnergyDensity" name="EDvoronoi" dynamic="e" static="ion0">
<spacegrid coord="voronoi"/>
</estimator>
The <reference_points/>
element provides a set of points for later
use in specifying the origin and coordinate axes needed to construct a
spatial histogramming grid. Several reference points on the surface of
the simulation cell (see Table 10), as well as the
positions of the ions (see the energydensity.static
attribute), are
made available by default. The reference points can be used, for
example, to construct a cylindrical grid along a bond with the origin on
the bond center.
reference_points
element:
parent elements:
estimator type=EnergyDensity
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
coord
\(^r\)text
Cartesian/cell
Specify coordinate system
body text: The body text is a line formatted list of points with labels
Additional information:
coord:
Ifcoord=cartesian
, labeled points are in Cartesian (x,y,z) format in units of Bohr. Ifcoord=cell
, then labeled points are in units of the simulation cell axes.body text:
The list of points provided in the body text are line formatted, with four entries per line (label coor1 coor2 coor3). A set of points referenced to the simulation cell is available by default (see Table 10). Ifenergydensity.static
is provided, the location of each individual ion is also available (e.g., ifenergydensity.static=ion0
, then the location of the first atom is available with label ion01, the second with ion02, etc.). All points can be used by label when constructing spatial histogramming grids (see the followingspacegrid
element) used to collect energy densities.





0 0 0 
Cell center 

\(a_1\) 
Cell axis 1 

\(a_2\) 
Cell axis 2 

\(a_3\) 
Cell axis 3 

\(a_1\)/2 
Cell face 1+ 

\(a_1\)/2 
Cell face 1 

\(a_2\)/2 
Cell face 2+ 

\(a_2\)/2 
Cell face 2 

\(a_3\)/2 
Cell face 3+ 

\(a_3\)/2 
Cell face 3 

\((a_1+a_2+a_3)/2\) 
Cell corner +,+,+ 

\((a_1+a_2a_3)/2\) 
Cell corner +,+, 

\((a_1a_2+a_3)/2\) 
Cell corner +,,+ 

\((a_1+a_2+a_3)/2\) 
Cell corner ,+,+ 

\((a_1a_2a_3)/2\) 
Cell corner +,, 

\((a_1+a_2a_3)/2\) 
Cell corner ,+, 

\((a_1a_2+a_3)/2\) 
Cell corner ,,+ 

\((a_1a_2a_3)/2\) 
Cell corner ,, 
Table 8 Reference points available by default. Vectors \(a_1\), \(a_2\), and \(a_3\) refer to the simulation cell axes. The representation of the cell is centered around zero
.
The <spacegrid/>
element is used to specify a spatial histogramming
grid for the energy density. Grids are constructed based on a set of,
potentially nonorthogonal, userprovided coordinate axes. The axes are
based on information available from reference_points
. Voronoi grids
are based only on nearest neighbor distances between electrons and ions.
Any number of space grids can be provided to a single energy density
estimator.
spacegrid
element:
parent elements:
estimator type=EnergyDensity
child elements:
origin, axis
attributes:
Name
Datatype
Values
Default
Description
coord
\(^r\)text
Cartesian
Specify coordinate system
cylindrical
spherical
Voronoi
The <origin/>
element gives the location of the origin for a
nonVoronoi grid.
Additional information:
p1/p2/fraction:
The location of the origin is set top1+fraction*(p2p1)
. If onlyp1
is provided, the origin is atp1
.
origin
element:
parent elements:
spacegrid
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
p1
\(^r\)text
reference_point.label
Select end point
p2
\(^o\)text
reference_point.label
Select end point
fraction
\(^o\)real
0
Interpolation fraction
The <axis/>
element represents a coordinate axis used to construct the, possibly curved, coordinate system for the histogramming grid. Three <axis/>
elements must be provided to a nonVoronoi <spacegrid/>
element.
axis
element:
parent elements:
spacegrid
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
label
\(^r\)text
See below
Axis/dimension label
grid
\(^r\)text
“0 1”
Grid ranges/intervals
p1
\(^r\)text
reference_point.label
Select end point
p2
\(^o\)text
reference_point.label
Select end point
scale
\(^o\)real
Interpolation fraction
Additional information:
label:
The allowed set of axis labels depends on the coordinate system (i.e.,spacegrid.coord
). Labels arex/y/z
forcoord=cartesian
,r/phi/z
forcoord=cylindrical
,r/phi/theta
forcoord=spherical
.p1/p2/scale:
The axis vector is set top1+scale*(p2p1)
. If onlyp1
is provided, the axis vector isp1
.grid:
The grid specifies the histogram grid along the direction specified bylabel
. The allowed grid points fall in the range [1,1] forlabel=x/y/z
or [0,1] forr/phi/theta
. A grid of 10 evenly spaced points between 0 and 1 can be requested equivalently bygrid="0 (0.1) 1"
orgrid="0 (10) 1."
Piecewise uniform grids covering portions of the range are supported, e.g.,grid="0.7 (10) 0.0 (20) 0.5."
Note that
grid
specifies the histogram grid along the (curved) coordinate given bylabel
. The axis specified byp1/p2/scale
does not correspond onetoone withlabel
unlesslabel=x/y/z
, but the full set of axes provided defines the (sheared) space on top of which the curved (e.g., spherical) coordinate system is built.
One body density matrix¶
The Nbody density matrix in DMC is \(\hat{\rho}_N=\left{\Psi_{T}}\rangle{}\langle{\Psi_{FN}}\right\) (for VMC, substitute \(\Psi_T\) for \(\Psi_{FN}\)). The one body reduced density matrix (1RDM) is obtained by tracing out all particle coordinates but one:
In this formula, the sum is over all electron indices and \(Tr_{R_n}(*)\equiv\int dR_n\langle{R_n}\left{*}\right{R_n}\rangle\) with \(R_n=[r_1,...,r_{n1},r_{n+1},...,r_N]\). When the sum is restricted over spinup or spindown electrons, one obtains a density matrix for each spin species. The 1RDM computed by is partitioned in this way.
In real space, the matrix elements of the 1RDM are
A more efficient and compact representation of the 1RDM is obtained by expanding in the SPOs obtained from a HartreeFock or DFT calculation, \(\{\phi_i\}\):
The integration over \(r'\) in (41) is inefficient when one is also interested in obtaining matrices involving energetic quantities, such as the energy density matrix of [KKR14] or the related (and more well known) generalized Fock matrix. For this reason, an approximation is introduced as follows:
For VMC, FNDMC, FPDMC, and RNDMC this formula represents an exact sampling of the 1RDM corresponding to \(\hat{\rho}_N^\dagger\) (see appendix A of [KKR14] for more detail).
estimtor type=dm1b
element:
parent elements:
hamiltonian, qmc
child elements:
None
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
dm1b
Must be dm1b
name
\(^r\)text
anything
Unique name for estimator
parameters:
Name
Datatype
Values
Default
Description
basis
\(^r\)text array
sposet.name(s)
Orbital basis
integrator
\(^o\)text
uniform_grid uniform density
uniform_grid
Integration method
evaluator
\(^o\)text
loop/matrix
loop
Evaluation method
scale
\(^o\)real
\(0<scale<1\)
1.0
Scale integration cell
center
\(^o\)real array(3)
any point
Center of cell
points
\(^o\)integer
\(>0\)
10
Grid points in each dim
samples
\(^o\)integer
\(>0\)
10
MC samples
warmup
\(^o\)integer
\(>0\)
30
MC warmup
timestep
\(^o\)real
\(>0\)
0.5
MC time step
use_drift
\(^o\)boolean
yes/no
no
Use drift in VMC
check_overlap
\(^o\)boolean
yes/no
no
Print overlap matrix
check_derivatives
\(^o\)boolean
yes/no
no
Check density derivatives
acceptance_ratio
\(^o\)boolean
yes/no
no
Print accept ratio
rstats
\(^o\)boolean
yes/no
no
Print spatial stats
normalized
\(^o\)boolean
yes/no
yes
basis
comes norm’ed
volume_normed
\(^o\)boolean
yes/no
yes
basis
norm is volume
energy_matrix
\(^o\)boolean
yes/no
no
Energy density matrix
Additional information:
name:
Density matrix results appear instat.h5
files labeled according toname
.basis:
Listsposet.name
’s. The total set of orbitals contained in allsposet
’s comprises the basis (subspace) onto which the one body density matrix is projected. This set of orbitals generally includes many virtual orbitals that are not occupied in a single reference Slater determinant.integrator:
Select the method used to perform the additional single particle integration. Options areuniform_grid
(uniform grid of points over the cell),uniform
(uniform random sampling over the cell), anddensity
(Metropolis sampling of approximate density, \(\sum_{b\in \texttt{basis}}\left{\phi_b}\right^2\), is not well tested, please check results carefully!). Depending on the integrator selected, different subsets of the other input parameters are active.evaluator:
Select forloop or matrix multiply implementations. Matrix is preferred for speed. Both implementations should give the same results, but please check as this has not been exhaustively tested.scale:
Resize the simulation cell by scale for use as an integration volume (active forintegrator=uniform/uniform_grid
).center:
Translate the integration volume to center at this point (active forintegrator=uniform/ uniform_grid
). Ifcenter
is not provided, the scaled simulation cell is used as is.points:
Number of grid points in each dimension forintegrator=uniform_grid
. For example,points=10
results in a uniform \(10 \times 10 \times 10\) grid over the cell.samples:
Sets the number of MC samples collected for each step (active forintegrator=uniform/ density
).warmup:
Number of warmup Metropolis steps at the start of the run before data collection (active forintegrator=density
).timestep:
Driftdiffusion time step used in Metropolis sampling (active forintegrator=density
).use_drift:
Enable drift in Metropolis sampling (active forintegrator=density
).check_overlap:
Print the overlap matrix (computed via simple Riemann sums) to the log, then abort. Note that subsequent analysis based on the 1RDM is simplest if the input orbitals are orthogonal.check_derivatives:
Print analytic and numerical derivatives of the approximate (sampled) density for several sample points, then abort.acceptance_ratio:
Print the acceptance ratio of the density sampling to the log for each step.rstats:
Print statistical information about the spatial motion of the sampled points to the log for each step.normalized:
Declare whether the inputted orbitals are normalized or not. Ifnormalized=no
, direct Riemann integration over a \(200 \times 200 \times 200\) grid will be used to compute the normalizations before use.volume_normed:
Declare whether the inputted orbitals are normalized to the cell volume (default) or not (a norm of 1.0 is assumed in this case). Currently, Bspline orbitals coming from QE and HEG planewave orbitals native to QMCPACK are known to be volume normalized.energy_matrix:
Accumulate the one body reduced energy density matrix, and write it tostat.h5
. This matrix is not covered in any detail here; the interested reader is referred to [KKR14].
<estimator type="dm1b" name="DensityMatrices">
<parameter name="basis" > spo_u spo_uv </parameter>
<parameter name="evaluator" > matrix </parameter>
<parameter name="integrator" > uniform_grid </parameter>
<parameter name="points" > 4 </parameter>
<parameter name="scale" > 1.0 </parameter>
<parameter name="center" > 0 0 0 </parameter>
</estimator>
<estimator type="dm1b" name="DensityMatrices">
<parameter name="basis" > spo_u spo_uv </parameter>
<parameter name="evaluator" > matrix </parameter>
<parameter name="integrator" > uniform </parameter>
<parameter name="samples" > 64 </parameter>
<parameter name="scale" > 1.0 </parameter>
<parameter name="center" > 0 0 0 </parameter>
</estimator>
<estimator type="dm1b" name="DensityMatrices">
<parameter name="basis" > spo_u spo_uv </parameter>
<parameter name="evaluator" > matrix </parameter>
<parameter name="integrator" > density </parameter>
<parameter name="samples" > 64 </parameter>
<parameter name="timestep" > 0.5 </parameter>
<parameter name="use_drift" > no </parameter>
</estimator>
<sposet_builder type="bspline" href="../dft/pwscf_output/pwscf.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" meshfactor="1.0" gpu="no" precision="single">
<sposet type="bspline" name="spo_u" group="0" size="4"/>
<sposet type="bspline" name="spo_d" group="0" size="2"/>
<sposet type="bspline" name="spo_uv" group="0" index_min="4" index_max="10"/>
</sposet_builder>
<sposet_builder type="bspline" href="../dft/pwscf_output/pwscf.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" meshfactor="1.0" gpu="no" precision="single">
<sposet type="bspline" name="spo_u" group="0" size="4"/>
<sposet type="bspline" name="spo_d" group="0" size="2"/>
<sposet type="bspline" name="dm_basis" size="50" spindataset="0"/>
</sposet_builder>
ForwardWalking Estimators¶
Forward walking is a method for sampling the pure fixednode distribution \(\langle \Phi_0  \Phi_0\rangle\). Specifically, one multiplies each walker’s DMC mixed estimate for the observable \(\mathcal{O}\), \(\frac{\mathcal{O}(\mathbf{R})\Psi_T(\mathbf{R})}{\Psi_T(\mathbf{R})}\), by the weighting factor \(\frac{\Phi_0(\mathbf{R})}{\Psi_T(\mathbf{R})}\). As it turns out, this weighting factor for any walker \(\mathbf{R}\) is proportional to the total number of descendants the walker will have after a sufficiently long projection time \(\beta\).
To forward walk on an observable, declare a generic forwardwalking
estimator within a <hamiltonian>
block, and then specify the
observables to forward walk on and the forwardwalking parameters. Here
is a summary.
estimator type=ForwardWalking
element:
parent elements:
hamiltonian, qmc
child elements:
Observable
attributes:
Name
Datatype
Values
Default
Description
type
\(^r\)text
ForwardWalking
Must be “ForwardWalking”
name
\(^r\)text
anything
any
Unique name for estimator
Observable
element:
parent elements:
estimator, hamiltonian, qmc
child elements:
None
Name
Datatype
Values
Default
Description
name
\(^r\)text
anything
any
Registered name of existing estimator on which to forward walk
max
\(^r\)integer
\(>0\)
Maximum projection time in steps (
max
\(=\beta/\tau\))
frequency
\(^r\)text
\(\geq 1\)
Dump data only for every
frequency
th toscalar.dat
file
Additional information:
Cost: Because histories of observables up to
max
time steps have to be stored, the memory cost of storing the nonforwardwalked observables variables should be multiplied by \(\texttt{max}\). Although this is not an issue for items such as potential energy, it could be prohibitive for observables such as density, forces, etc.Naming Convention: Forwardwalked observables are automatically named
FWE_name_i
, wherei
is the forwardwalked expectation value at time stepi
, andname
is whatever name appears in the<Observable>
block. This is also how it will appear in thescalar.dat
file.
In the following example case, QMCPACK forward walks on the potential energy for 300 time steps and dumps the forwardwalked value at every time step.
<estimator name="fw" type="ForwardWalking">
<Observable name="LocalPotential" max="300" frequency="1"/>
<! Additional Observable blocks go here >
</estimator>
ChiesaCeperleyZhang Force Estimators¶
All force estimators implemented in QMCPACK are invoked with type="Force"
.
The mode
determines the specific estimator to be used. Currently,
QMCPACK supports ChiesaCeperleyZhang (CCZ) allelectron and
AssarafCaffarel ZeroVariance ZeroBias (AC) force estimators. We’ll discuss
the CCZ estimator in this section, and the AC estimator in the following section.
Without loss of generality, the CCZ estimator for the zcomponent of the force on an ion centered at the origin is given by the following expression:
Z is the ionic charge, \(M\) is the degree of the smoothing polynomial, \(\mathcal{R}\) is a realspace cutoff of the sphere within which the bareforce estimator is smoothed, and \(c_\ell\) are predetermined coefficients. These coefficients are chosen to minimize the weighted mean square error between the bare force estimate and the swave filtered estimator. Specifically,
Here, \(m\) is the weighting exponent, \(f_z(r)\) is the unfiltered radial force density for the z force component, and \(\tilde{f}_z(r)\) is the smoothed polynomial function for the same force density.
Currently, open and periodic boundary conditions are supported but for allelectron calculations only.
The reader is invited to refer to the original paper for a more thorough explanation of the methodology, but with the notation in hand, QMCPACK takes the following parameters.
estimator type=Force
element:
parent elements:
hamiltonian, qmc
child elements:
parameter
attributes:
Name
Datatype
Values
Default
Description
mode
\(^o\)text
See above
bare
Select estimator type
lrmethod
\(^o\)text
ewald or srcoul
ewald
Select longrange potential breakup method
type
\(^r\)text
Force
Must be “Force”
name
\(^o\)text
Anything
ForceBase
Unique name for this estimator
pbc
\(^o\)boolean
yes/no
yes
Using periodic BCs or not
addionion
\(^o\)boolean
yes/no
no
Add the ionion force contribution to output force estimate
parameters:
Name
Datatype
Values
Default
Description
rcut
\(^o\)real
\(>0\)
1.0
Realspace cutoff \(\mathcal{R}\) in bohr
nbasis
\(^o\)integer
\(>0\)
2
Degree of smoothing polynomial \(M\)
weightexp
\(^o\)integer
\(>0\)
2
\(\chi^2\) weighting exponent :math`m`
Additional information:
Naming Convention: The unique identifier
name
is appended withname_X_Y
in thescalar.dat
file, whereX
is the ion ID number andY
is the component ID (an integer with x=0, y=1, z=2). All force components for all ions are computed and dumped to thescalar.dat
file.Longrange breakup: With periodic boundary conditions, it is important to converge the lattice sum when calculating Coulomb contribution to the forces. As a quick test, increase the
LR_dim_cutoff
parameter until ionion forces are converged. The Ewald method converges more slowly than optimized method, but the optimized method can break down in edge cases, eg. too largeLR_dim_cutoff
.Miscellaneous: Usually, the default choice of
weightexp
is sufficient. Different combinations ofrcut
andnbasis
should be tested though to minimize variance and bias. There is, of course, a tradeoff, with largernbasis
and smallerrcut
leading to smaller biases and larger variances.
The following is an example use case.
<simulationcell>
...
<parameter name="LR_handler"> opt_breakup_original </parameter>
<parameter name="LR_dim_cutoff"> 20 </parameter>
</simulationcell>
<hamiltonian>
<estimator name="F" type="Force" mode="cep" addionion="yes">
<parameter name="rcut">0.1</parameter>
<parameter name="nbasis">4</parameter>
<parameter name="weightexp">2</parameter>
</estimator>
</hamiltonian>
AssarafCaffarel Force Estimators¶
*WARNING: The following estimator formally has infinite variance. You MUST do something to mitigate this in postprocessing or during the run before publishing.*
QMCPACK has an implementation of force estimation using the AssarafCaffarel ZeroVariance ZeroBias method [TCK21]. This has the desirable property that as the trial wave function and trial wave function derivative become exact, the estimator itself becomes an exact estimate of the force and the variance of the estimator goes to ero–much like the local energy. In practice, the estimator is usually significantly more accurate and has much lower variance than the bare HellmanFeynman estimator.
Currently, this is the only recommended way to estimate forces for systems with nonlocal pseudopotentials.
The zerovariance, zerobias force estimator is given by the following expression:
The first term is the zerovariance force estimator, given by the following.
The first term is the bare “HellmanFeynman” term (denoted “hf” in output), and the second is a fluctuation cancelling zerovariance term (called “pulay” in output). This splitting allows the user to investigate the individual contributions to the force estimator, but we recommend always adding “hf” and “pulay” terms unless there is a compelling reason to do otherwise.
The second term is the “zerobias” term:
Because knowledge of \(\langle E_L \rangle\) is needed to compute the zerobias term, QMCPACK returns \(E_L(\mathbf{R}) \ln \Psi_T\) (called “Ewfngrad” in output), and \(\ln \Psi_T\) (called “wfngrad” in output), which in conjunction with the local energy, allows one to construct the zerobias term in postprocessing.
There is an initial implementation of spacewarp variance reduction that is accessible to the enduser, subject to the caveat that evaluation of these terms is currently slow (derivatives of local energy are computed with finite differences, rather than analytically). If the spacewarp option is enabled, the following term is added to the zerovariance force estimator:
The variance reduction with this term is observed to be rather large. A faster, more efficient implementation of this term will be available in a future QMCPACK release.
The following term is added to the wave function gradient:
Currently, there is only one choice for damping function \(\omega_I(\mathbf{r})\). This is given by:
We use \(F(r)=r^{4}\) for the real space damping.
Finally, the estimator provides two methods to evaluate the necessary derivatives of the wave function and Hamiltonian. The first is a straightforward analytic differentiation of all required terms. While mathematically transparent, this algorithm has poor scaling with system size. The second utilizes the fastderivative algorithm of Assaraf, Moroni, and Filippi [FAM16], which has a smaller computational prefactor and at least an O(N) speedup over the legacy implementation. Both of these methods are accessible with appropraite flags.
estimator type=Force
element:
parent elements:
hamiltonian, qmc
child elements:
none
attributes:
Name
Datatype
Values
Default
Description
mode
\(^o\)text
acforce
Required to use ACForce estimator
type
\(^r\)text
Force
Must be “Force”
name
\(^o\)text
Anything
ForceBase
Unique name for this estimator
spacewarp
\(^o\)text
yes/no
no
Add spacewarp variance reduction terms
fast_derivatives
\(^o\)text
yes/no
no
Use Filippi fast derivative algorithm
Additional information:
Naming Convention: The unique identifier
name
is appended with a term label (hf
,pulay
,Ewfngrad
, orwfgrad
)name_term_X_Y
in thescalar.dat
file, whereX
is the ion ID number andY
is the component ID (an integer with x=0, y=1, z=2). All force components for all ions are computed and dumped to thescalar.dat
file.Note: The fast force algorithm returns the total derivative of the local energy, and does not make the split between “HellmanFeynman” and zerovariance terms like the legacy force implementation does. As such, expect
name_pulay_X_Y
to be zero iffast_derivatives="yes"
. However, this will be identical to the sum of HellmanFeynman and zerovariance terms in the legacy implementation.
The following is an example use case.
<hamiltonian>
<estimator name="F" type="Force" mode="acforce" fast_derivatives="yes" spacewarp="no"/>
</hamiltonian>
Stress estimators¶
QMCPACK takes the following parameters.
parent elements:
hamiltonian
attributes:
Name
Datatype
Values
Default
Description
mode
\(^r\)text
stress
bare
Must be “stress”
type
\(^r\)text
Force
Must be “Force”
source
\(^r\)text
ion0
Name of ion particleset
name
\(^o\)text
Anything
ForceBase
Unique name for this estimator
addionion
\(^o\)boolean
yes/no
no
Add the ionion stress contribution to output
Additional information:
Naming Convention: The unique identifier
name
is appended withname_X_Y
in thescalar.dat
file, whereX
andY
are the component IDs (an integer with x=0, y=1, z=2).Longrange breakup: With periodic boundary conditions, it is important to converge the lattice sum when calculating Coulomb contribution to the forces. As a quick test, increase the
LR_dim_cutoff
parameter until ionion stresses are converged. Check using QE “Ewald contribution”, for example. The stress estimator is implemented only with the Ewald method.
The following is an example use case.
<simulationcell>
...
<parameter name="LR_handler"> ewald </parameter>
<parameter name="LR_dim_cutoff"> 45 </parameter>
</simulationcell>
<hamiltonian>
<estimator name="S" type="Force" mode="stress" source="ion0"/>
</hamiltonian>
 CCMH06
Simone Chiesa, David M. Ceperley, Richard M. Martin, and Markus Holzmann. Finitesize error in manybody simulations with longrange interactions. Phys. Rev. Lett., 97:076404, August 2006. doi:10.1103/PhysRevLett.97.076404.
 FAM16
Claudia Filippi, Roland Assaraf, and Saverio Moroni. Simple formalism for efficient derivatives and multideterminant expansions in quantum monte carlo. The Journal of Chemical Physics, 144(19):194105, 2016. URL: https://doi.org/10.1063/1.4948778, arXiv:https://doi.org/10.1063/1.4948778, doi:10.1063/1.4948778.
 KKR14(1,2,3)
Jaron T. Krogel, Jeongnim Kim, and Fernando A. Reboredo. Energy density matrix formalism for interacting quantum systems: quantum monte carlo study. Phys. Rev. B, 90:035125, July 2014. doi:10.1103/PhysRevB.90.035125.
 KYKC13
Jaron T. Krogel, Min Yu, Jeongnim Kim, and David M. Ceperley. Quantum energy density: improved efficiency for quantum monte carlo calculations. Phys. Rev. B, 88:035137, July 2013. doi:10.1103/PhysRevB.88.035137.
 MSC91
Lubos Mitas, Eric L. Shirley, and David M. Ceperley. Nonlocal pseudopotentials and diffusion monte carlo. The Journal of Chemical Physics, 95(5):3467–3475, 1991. doi:10.1063/1.460849.
 NC95
Vincent Natoli and David M. Ceperley. An optimized method for treating longrange potentials. Journal of Computational Physics, 117(1):171–178, 1995. URL: http://www.sciencedirect.com/science/article/pii/S0021999185710546, doi:10.1006/jcph.1995.1054.
 TCK21
Juha Tiihonen, Raymond C. Clay, and Jaron T. Krogel. Toward quantum monte carlo forces on heavier ions: scaling properties. The Journal of Chemical Physics, 154(20):204111, 2021. URL: https://doi.org/10.1063/5.0052266, arXiv:https://doi.org/10.1063/5.0052266, doi:10.1063/5.0052266.
 ZBMAlfe19
Andrea Zen, Jan Gerit Brandenburg, Angelos Michaelides, and Dario Alfè. A new scheme for fixed node diffusion quantum monte carlo with pseudopotentials: improving reproducibility and reducing the trialwavefunction bias. The Journal of Chemical Physics, 151(13):134105, October 2019. doi:10.1063/1.5119729.