AuxiliaryField Quantum Monte Carlo
The AFQMC method is an orbitalspace formulation of the imaginarytime propagation algorithm. We refer the reader to one of the review articles on the method [PZ04, Zha13, ZK03] for a detailed description of the algorithm. It uses the HubbardStratonovich transformation to express the imaginarytime propagator, which is inherently a 2body operator, as an integral over 1body propagators, which can be efficiently applied to an arbitrary Slater determinant. This transformation allows us to represent the interacting manybody system as an average over a noninteracting system (e.g., Slater determinants) in a timedependent fluctuating external field (the Auxiliary fields). The walkers in this case represent nonorthogonal Slater determinants, whose time average represents the desired quantum state. QMCPACK currently implements the phaseless AFQMC algorithm of Zhang and Krakauer [ZK03], where a trial wavefunction is used to project the simulation to the real axis, controlling the fermionic sign problem at the expense of a bias. This approximation is similar in spirit to the fixednode approximation in realspace DMC but applied in the Hilbert space where the AFQMC random walk occurs.
Input
The input for an AFQMC calculation is fundamentally different to the
input for other realspace algorithms in QMCPACK. The main source of
input comes from the Hamiltonian matrix elements in an appropriate
single particle basis. This must be evaluated by an external code and
saved in a format that QMCPACK can read. More details about file formats
follow. The input file has six basic xmlblocks: AFQMCInfo
,
Hamiltonian
, Wavefunction
, WalkerSet
, Propagator
, and
execute
. The first five define input structures required for various
types of calculations. The execute
block represents actual
calculations and takes as input the other blocks. Nonexecution blocks
are parsed first, followed by a second pass where execution blocks are
parsed (and executed) in order. Listing 51 shows an example of a
minimal input file for an AFQMC calculation.
table13
shows a brief description of the most
important parameters in the calculation. All xml sections contain a
“name” argument used to identify the resulting object within QMCPACK.
For example, in the example, multiple Hamiltonian objects with different
names can be defined. The one actually used in the calculation is the
one passed to “execute” as ham.
<?xml version="1.0"?>
<simulation method="afqmc">
<project id="Carbon" series="0"/>
<AFQMCInfo name="info0">
<parameter name="NMO">32</parameter>
<parameter name="NAEA">16</parameter>
<parameter name="NAEB">16</parameter>
</AFQMCInfo>
<Hamiltonian name="ham0" info="info0">
<parameter name="filename">fcidump.h5</parameter>
</Hamiltonian>
<Wavefunction name="wfn0" type="MSD" info="info0">
<parameter name="filetype">hdf5</parameter>
<parameter name="filename">wfn.h5</parameter>
</Wavefunction>
<WalkerSet name="wset0">
<parameter name="walker_type">closed</parameter>
</WalkerSet>
<Propagator name="prop0" info="info0">
</Propagator>
<execute wset="wset0" ham="ham0" wfn="wfn0" prop="prop0" info="info0">
<parameter name="timestep">0.005</parameter>
<parameter name="blocks">10000</parameter>
<parameter name="nWalkers">20</parameter>
<Estimator name="back_propagation">
<parameter name="naverages">4</parameter>
<parameter name="nsteps">400</parameter>
<parameter name="path_restoration">true</parameter>
<onerdm/>
<diag2rdm/>
<twordm/>
<ontop2rdm/>
<realspace_correlators/>
<correlators/>
<genfock/>
</Estimator>
</execute>
</simulation>
The following list includes all input sections for AFQMC calculations, along with a detailed explanation of accepted parameters. Since the code is under active development, the list of parameters and their interpretation might change in the future.
AFQMCInfo
: Input block that defines basic information about the
calculation. It is passed to all other input blocks to propagate the
basic information: <AFQMCInfo name="info0">
NMO. Number of molecular orbitals, i.e., number of states in the single particle basis.
NAEA. Number of active electronsalpha, i.e., number of spinup electrons.
NAEB. Number of active electronsbeta, i.e., number of spindown electrons.
Hamiltonian
: Controls the object that reads, stores, and manages the
hamiltonian
.
<Hamiltonian name="ham0" type="SparseGeneral" info="info0">
filename. Name of file with the
Hamiltonian
. This is a required parameter.cutoff_1bar. Cutoff applied to integrals during reading. Any term in the Hamiltonian smaller than this value is set to zero. (For filetype=“hdf5”, the cutoff is applied only to the 2electron integrals). Default: 1e8
cutoff_decomposition. Cutoff used to stop the iterative cycle in the generation of the Cholesky decomposition of the 2electron integrals. The generation of Cholesky vectors is stopped when the maximum error in the diagonal reaches this value. In case of an eigenvalue factorization, this becomes the cutoff applied to the eigenvalues. Only eigenvalues above this value are kept. Default: 1e6
nblocks. This parameter controls the distribution of the 2electron integrals among processors. In the default behavior (nblocks=1), all nodes contain the entire list of integrals. If nblocks \(>\) 1, the of nodes in the calculation will be split in nblocks groups. Each node in a given group contains the same subset of integrals and subsequently operates on this subset during any further operation that requires the hamiltonian. The maximum number of groups is NMO. Currently only works for filetype=“hdf5” and the file must contain integrals. Not yet implemented for input hamiltonians in the form of Cholesky vectors or for ASCII input. Coming soon! Default: No distribution
printEig. If “yes”, prints additional information during the Cholesky decomposition. Default: no
fix_2eint. If this is set to “yes”, orbital pairs that are found not to be positive definite are ignored in the generation of the Cholesky factorization. This is necessary if the 2electron integrals are not positive definite because of roundoff errors in their generation. Default: no
Wavefunction
: controls the object that manages the trial
wavefunctions. This block expects a list of xmlblocks defining actual
trial wavefunctions for various roles.
<Wavefunction name="wfn0" type="MSD/PHMSD" info="info0">
filename. Name of file with wavefunction information.
cutoff. cutoff applied to the terms in the calculation of the local energy. Only terms in the Hamiltonian above this cutoff are included in the evaluation of the energy. Default: 1e6
nnodes. Defines the parallelization of the local energy evaluation and the distribution of the
Hamiltonian
matrix (not to GPU)nbatch_qr. This turns on(>=1)/off(==0) batched QR calculation. 1 means all the walkers in the batch. Default: 0 (CPU) / 1 (GPU)
WalkerSet
: Controls the object that handles the set of walkers.
<WalkerSet name="wset0">
walker_type. Type of walker set: closed or collinear. Default: collinear
pop_control. Population control algorithm. Options: “simple”: Uses a simple branching scheme with a fluctuating population. Walkers with weight above max_weight are split into multiple walkers of weight reset_weight. Walkers with weight below min_weight are killed with probability (weight/min_weight); “pair”: Fixedpopulation branching algorithm, based on QWalk’s branching algorithm. Pairs of walkers with weight above/below max_weight/min_weight are combined into 2 walkers with weights equal to \((w_1+w_2)/2\). The probability of replicating walker w1 (larger weight) occurs with probability \(w_1/(w_1+w_2)\), otherwise walker w2 (lower weight) is replicated; “comb”: Fixedpopulation branching algorithm based on the Comb method. Will be available in the next release. Default: “pair”
min_weight. Weight at which walkers are possibly killed (with probability weight/min_weight). Default: 0.05
max_weight. Weight at which walkers are replicated. Default: 4.0
reset_weight. Weight to which replicated walkers are reset to. Default: 1.0
Propagator
: Controls the object that manages the propagators.
<Propagator name="prop0" info="info0">
cutoff. Cutoff applied to Cholesky vectors. Elements of the Cholesky vectors below this value are set to zero. Only meaningful with sparse hamiltonians. Default: 1e6
substractMF. If “yes”, apply meanfield subtraction based on the ImpSamp trial wavefunction. Must set to “no” to turn it off. Default: yes
vbias_bound. Upper bound applied to the vias potential. Components of the vias potential above this value are truncated there. The bound is currently applied to \(\sqrt{\tau} v_{bias}\), so a larger value must be used as either the time step or the fluctuations increase (e.g. from running a larger system or using a poor trial wavefunction). Default: 3.0
apply_constrain. If “yes”, apply the phaseless constrain to the walker propagation. Currently, setting this to “no” produces unknown behavior, since free propagation algorithm has not been tested. Default: yes
hybrid. If “yes”, use hybrid propagation algorithm. This propagation scheme doesn’t use the local energy during propagation, leading to significant speed ups when its evaluation cost is high. The local energy of the ImpSamp trial wavefunction is never evaluated. To obtain energy estimates in this case, you must define an Estimator xmlblock with the
Wavefunction
block. The local energy of this trial wavefunction is evaluated and printed. It is possible to use a previously defined trial wavefunction in the Estimator block, just set its “name” argument to the name of a previously defined wavefunction. In this case, the same object is used for both roles. Default: nonnodes. Controls the parallel propagation algorithm. If nnodes \(>\) 1, the nodes in the simulation are split into groups of nnodes nodes, each group working collectively to propagate their walkers. Default: 1 (Serial algorithm)
nbatch. This turns on(>=1)/off(==0) batched calculation of density matrices and overlaps. 1 means all the walkers in the batch. Default: 0 (CPU) / 1 (GPU)
nbatch_qr. This turns on(>=1)/off(==0) batched QR calculation. 1 means all the walkers in the batch. Default: 0 (CPU) / 1 (GPU)
execute
: Defines an execution region.
<execute wset="wset0" ham="ham0" wfn="wfn0" prop="prop0" info="info0">
nWalkers. Initial number of walkers per core group (see ncores). This sets the number of walkers for a given group of “ncores” on a node; the total number of walkers in the simulation depends on the total number of nodes and on the total number of cores on a node in the following way: \(\#_walkers_total = nWalkers * \#_nodes * \#_cores_total / ncores\). Default: 5
timestep. Time step in 1/a.u. Default: 0.01
blocks. Number of blocks. Slow operations occur once per block (e.g., write to file, slow observables, checkpoints), Default: 100
step. Number of steps within a block. Operations that occur at the step level include load balance, orthogonalization, branching, etc. Default: 1
substep. Number of substeps within a step. Only walker propagation occurs in a substep. Default: 1
ortho. Number of steps between orthogonalization. Default: 1
ncores. Number of nodes in a task group. This number defines the number of cores on a node that share the parallel work associated with a distributed task. This number is used in the
Wavefunction
andPropagator
task groups. The walker sets are shares by the ncores on a given node in the task group.checkpoint. Number of blocks between checkpoint files are generated. If a value smaller than 1 is given, no file is generated. If hdf_write_file is not set, a default name is used. Default: 0
hdf_write_file. If set (and checkpoint>0), a checkpoint file with this name will be written.
hdf_read_file. If set, the simulation will be restarted from the given file.
Within the Estimators
xml block has an argument name: the type
of estimator we want to measure. Currently available estimators include:
“basic”, “energy”, “mixed_one_rdm”, and “back_propagation”.
The basic estimator has the following optional parameters:
timers. print timing information. Default: true
The back_propagation estimator has the following parameters:
ortho. Number of backpropagation steps between orthogonalization. Default: 10
nsteps. Maximum number of backpropagation steps. Default: 10
naverages. Number of back propagation calculations to perform. The number of steps will be chosed equally distributed in the range 0,nsteps. Default: 1
block_size. Number of blocks to use in the internal average of the back propagated estimator. This is used to block data and reduce the size of the output. Default: 1
nskip. Number of blocks to skip at the start of the calculation for equilibration purposes. Default: 0
path_restoration. Use full path restoration. Can result in better back propagated results. Default false.
The following observables can be computed with the back_propagated estimator
onerdm. Oneparticle reduced density matrix.
twordm. Full Twoparticle reduced density matrix.
diag2rdm. Diagonal part of the twoparticle reduced density matrix.
ontop2rdm. On top twoparticle reduced density matrix.
realspace_correlators. ChargeCharge, and spinspin correlation functions in real space.
correlators. ChargeCharge, and spinspin correlation functions in real space centered about atomic sites.
genfock. Generalized Fock matrix.
Real space correlation functions require a real space grid. Details coming soon..
Hamiltonian File formats
QMCPACK offers three factorization approaches which are appropriate in different settings. The most generic approach implemented is based on the modifiedCholesky factorization [ADVFerre+09, BL77, KdMerasP03, PKVZ11, PZK13] of the ERI tensor:
where the sum is truncated at \(N_{\mathrm{chol}} = x_c M\), \(x_c\) is typically between \(5\) and \(10\), \(M\) is the number of basis functions and we have assumed that the singleparticle orbitals are in general complex. The storage requirement is thus naively \(\mathcal{O}(M^3)\). Note we follow the usual definition of \(v_{pqrs} = \langle pq  rs \rangle = (prqs)\). With this form of factorization QMCPACK allows for the integrals to be stored in either dense or sparse format.
The dense case is the simplest and is only implemented for Hamiltonians with real integrals (and basis functions, i.e. not the homegeneous electron gas which has complex orbitals but real integrals). The file format is given as follows:
$ h5dump n afqmc.h5
HDF5 "afqmc.h5" {
FILE_CONTENTS {
group /
group /Hamiltonian
group /Hamiltonian/DenseFactorized
dataset /Hamiltonian/DenseFactorized/L
dataset /Hamiltonian/dims
dataset /Hamiltonian/hcore
dataset /Hamiltonian/Energies
}
}
where the datasets are given by the following
/Hamiltonian/DenseFactorized/L
Contains the \([M^2,N_\mathrm{nchol}]\) dimensional matrix representatation of \(L_{pr,n}\)./Hamiltonian/dims
Descriptor array of length 8 containing \([0,0,0,M,N_\alpha,N_\beta,0,N_\mathrm{nchol}]\). Note that \(N_\alpha\) and \(N_\beta\) are somewhat redundant and will be read from the input file and wavefunction. This allows for the Hamiltonian to be used with different (potentially spin polarized) wavefunctions./Hamiltonian/hcore
Contains the \([M,M]\) dimensional onebody Hamiltonian matrix elements \(h_{pq}\)./Hamiltonian/Energies
Array containing \([E_{II}, E_{\mathrm{core}}]\). \(E_{II}\) should contain ionion repulsion energy and any additional constant terms which have to be added to the total energy. \(E_{\mathrm{core}}\) is deprecated and not used.
Typically the Cholesky matrix is sparse, particularly if written in the nonorthogonal AO basis (not currently supported in QMCPACK). In this case only a small number of nonzero elements (denoted \(nnz\) below) need to be stored which can reduce the memory overhead considerably. Internally QMCPACK stores this matrix in the CSR format, and the HDF5 file format is reflective of this. For large systems and, more generally when running in parallel, it is convenient to chunk the writing/reading of the Cholesky matrix into blocks of size \([M^2,\frac{N_{\mathrm{chol}}}{N_{\mathrm{blocks}}}]\) (if interpreted as a dense array). This is achieved by writing these blocks to different data sets in the file. For the sparse case the Hamtiltonian file format is given as follows:
$ h5dump n afqmc.h5
HDF5 "afqmc.h5" {
FILE_CONTENTS {
group /
group /Hamiltonian
group /Hamiltonian/Factorized
dataset /Hamiltonian/Factorized/block_sizes
dataset /Hamiltonian/Factorized/index_0
dataset /Hamiltonian/Factorized/vals_0
dataset /Hamiltonian/ComplexIntegrals
dataset /Hamiltonian/dims
dataset /Hamiltonian/hcore
dataset /Hamiltonian/Energies
}
}
/Hamiltonian/Factorized/block_sizes
Contains the number of elements in each block of the sparse representation of the Cholesky matrix \(L_{pr,n}\). In this case there is 1 block./Hamiltonian/Factorized/index_0
\([2\times nnz]\) dimensional array, containing the indices of the nonzero values of \(L_{ik,n}\). The row indices are stored in the even entries, and the column indices in the odd entries./Hamiltonian/Factorized/vals_0
\([nnz]\) length array containing nonzero values of \(L_{pr,n}\) for chunk 0./Hamiltonian/dims
Descriptor array of length 8 containing \([0,nnz,N_{\mathrm{block}},M,N_\alpha,N_\beta,0,N_\mathrm{nchol}]\)./Hamiltonian/ComplexIntegrals
Length 1 array that specifies if integrals are complex valued. 1 for complex integrals, 0 for real integrals./Hamiltonian/hcore
Contains the \([M,M]\) dimensional onebody Hamiltonian matrix elements \(h_{pq}\). Due to its small size this is written as a dense 2Darray./Hamiltonian/Energies
Array containing \([E_{II}, E_{\mathrm{core}}]\). \(E_{II}\) should contain ionion repulsion energy and any additional constant terms which have to be added to the total energy. \(E_{\mathrm{core}}\) is deprecated and not used.
To reduce the memory overhead of storing the threeindex tensor we recently adapted the tensorhypercontraction [HPMartinez12, HPSMartinez12, PHMartinezS12] (THC) approach for use in AFQMCcite{MaloneISDF2019}. Within the THC approach we can approximate the orbital products entering the ERIs as
where \(\varphi_p(\mathbf{r})\) are the oneelectron orbitals and \(\mathbf{r}_\mu\) are a set of specially selected interpolating points, \(\zeta_\mu(\mathbf{r})\) are a set of interpolating vectors and \(N_\mu = x_\mu M\). We can then write the ERI tensor as a product of rank2 tensors
where
We also require the halfrotated versions of these quantities which live on a different set of \(\tilde{N}_\mu\) interpolating points \(\tilde{\mathbf{r}}_\mu\) (see [MZM19]). The file format for THC factorization is as follows:
$ h5dump n afqmc.h5
HDF5 "afqmc.h5" {
FILE_CONTENTS {
group /
group /Hamiltonian
group /Hamiltonian/THC
dataset /Hamiltonian/THC/Luv
dataset /Hamiltonian/THC/Orbitals
dataset /Hamiltonian/THC/HalfTransformedMuv
dataset /Hamiltonian/THC/HalfTransformedFullOrbitals
dataset /Hamiltonian/THC/HalfTransformedOccOrbitals
dataset /Hamiltonian/THC/dims
dataset /Hamiltonian/ComplexIntegrals
dataset /Hamiltonian/dims
dataset /Hamiltonian/hcore
dataset /Hamiltonian/Energies
}
}
/Hamiltonian/THC/Luv
Cholesky factorization of the \(M_{\mu\nu}\) matrix given in (74)./Hamiltonian/THC/Orbitals
\([M,N_\mu]\) dimensional array of orbitals evaluated at chosen interpolating points \(\varphi_i(\mathbf{r}_\mu)\)./Hamiltonian/THC/HalfTransformedMuv
\([\tilde{N}_\mu,\tilde{N}_\mu]\) dimensional array containing halftransformed \(\tilde{M}_{\mu\nu}\)./Hamiltonian/THC/HalfTransformedFullOrbitals
\([M,\tilde{N}_\mu]\) dimensional array containing orbital set computed at halftransformed interpolating points \(\varphi_i(\tilde{\mathbf{r}}_\mu)\)./Hamiltonian/THC/HalfTransformedOccOrbitals
\([N_\alpha+N_\beta,\tilde{N}_\mu]\) dimensional array containing halfrotated orbital set computed at halftransformed interpolating points \(\varphi_a(\tilde{\mathbf{r}}_\mu) = \sum_{p} A_{pa}^* \varphi_{p}(\tilde{\mathbf{r}}_\mu)\), where \(\mathbf{A}\) is the SlaterMatrix of the (currently singledeterminant) trial wavefunction./Hamiltonian/THC/dims
Descriptor array containing \([M, N_\mu, \tilde{N}_\mu]\)./Hamiltonian/ComplexIntegrals
Length 1 array that specifies if integrals are complex valued. 1 for complex integrals, 0 for real integrals./Hamiltonian/dims
Descriptor array of length 8 containing \([0,0,0,M,N_\alpha,N_\beta,0,0]\)./Hamiltonian/hcore
Contains the \([M,M]\) dimensional onebody Hamiltonian matrix elements \(h_{ij}\)./Hamiltonian/Energies
Array containing \([E_{II}, E_{\mathrm{core}}]\). \(E_{II}\) should contain ionion repulsion energy and any additional constant terms which have to be added to the total energy (such as the electronelectron interaction Madelung contribution of \(\frac{1}{2} N \xi )\). \(E_{\mathrm{core}}\) is deprecated and not used.
Finally, we have implemented an explicitly \(k\)point dependent factorization for periodic systems [MZM20, MZC19]
where \(\textbf{k}\), \(\textbf{k}'\) and \(\textbf{Q}\) are vectors in the first Brillouin zone. The onebody Hamiltonian is block diagonal in \(\textbf{k}\) and in (75) we have used momentum conservation \((\textbf{k}_p  \textbf{k}_r + \textbf{k}_q  \textbf{k}_s) = \textbf{G}\) with \(\textbf{G}\) being some vector in the reciprocal lattice of the simulation cell. The convention for the Cholesky matrix \(L^{\textbf{Q},\textbf{k}}_{pr,\gamma}\) is as follows: \(\textbf{k}_r = \textbf{k}_p  \textbf{Q}\), so the vector \(\textbf{k}\) labels the kpoint of the first band index, \(\textit{p}\), while the kpoint vector of the second band index, \(\textit{r}\), is given by \(\textbf{k}  \textbf{Q}\). Electron repulsion integrals at different \(\textbf{Q}\) vectors are zero by symmetry, resulting in a reduction in the number of required \(\mathbf{Q}\) vectors. For certain \(\textbf{Q}\) vectors that satisfy \(\textbf{Q} \ne \textbf{Q}\) (this is not satisfied at the origin and at high symmetry points on the edge of the 1BZ), we have \({L^{\textbf{Q},\textbf{k}}_{sq,\gamma}}^{*} = {L^{\textbf{Q},\textbf{k}\textbf{Q}}_{qs,\gamma}}\), which requires us to store Cholesky vectors for either one of the \((\textbf{Q},\textbf{Q})\) pair, but not both.
In what follows let \(m_{\mathbf{k}}\) denote the number of basis functions for basis functions of a given \(k\)point (these can in principle differ for different \(k\)points due to linear dependencies), \(n^{\alpha}_{\mathbf{k}}\) the number of \(\alpha\) electrons in a given \(k\)point and \(n_{\mathrm{chol}}^{\mathbf{Q}_n}\) the number of Cholesky vectors for momentum transfer \(\mathbf{Q}_n\). The file format for this factorization is as follows (for a \(2\times2\times2\) \(k\)point mesh, for denser meshes generally there will be far fewer symmetry inequivalent momentum transfer vectors than there are \(k\)points):
$ h5dump n afqmc.h5
HDF5 "afqmc.h5" {
FILE_CONTENTS {
group /
group /Hamiltonian
group /Hamiltonian/KPFactorized
dataset /Hamiltonian/KPFactorized/L0
dataset /Hamiltonian/KPFactorized/L1
dataset /Hamiltonian/KPFactorized/L2
dataset /Hamiltonian/KPFactorized/L3
dataset /Hamiltonian/KPFactorized/L4
dataset /Hamiltonian/KPFactorized/L5
dataset /Hamiltonian/KPFactorized/L6
dataset /Hamiltonian/KPFactorized/L7
dataset /Hamiltonian/NCholPerKP
dataset /Hamiltonian/MinusK
dataset /Hamiltonian/NMOPerKP
dataset /Hamiltonian/QKTok2
dataset /Hamiltonian/H1_kp0
dataset /Hamiltonian/H1_kp1
dataset /Hamiltonian/H1_kp2
dataset /Hamiltonian/H1_kp3
dataset /Hamiltonian/H1_kp4
dataset /Hamiltonian/H1_kp5
dataset /Hamiltonian/H1_kp6
dataset /Hamiltonian/H1_kp7
dataset /Hamiltonian/ComplexIntegrals
dataset /Hamiltonian/KPoints
dataset /Hamiltonian/dims
dataset /Hamiltonian/Energies
}
}
/Hamiltonian/KPFactorized/L[n]
This series of datasets store elements of the Cholesky tensors \(L[\mathbf{Q}_n,\mathbf{k},pr,n]\). Each data set is of dimension \([N_k,m_{\mathbf{k}}\times m_{\mathbf{k}'},n^{\mathbf{Q}_n}_\mathrm{chol}]\), where, again, \(k\) is the \(k\)point associated with basis function \(p\), the \(k\)point of basis function \(r\) is defined via the mappingQKtok2
./Hamiltonian/NCholPerKP
\(N_k\) length array giving number of Cholesky vectors per \(k\)point./Hamiltonian/MinusK
: \(N_k\) length array mapping a \(k\)point to its inverse: \(\mathbf{k}_i+\)MinusK[i]
\(= \mathbf{0} \mod \mathbf{G}\)./Hamiltonian/NMOPerKP
: \(N_k\) length array listing number of basis functions per \(k\)point./Hamiltonian/QKTok2
: \([N_k,N_k]\) dimensional array.QKtok2[i,j]
yields the \(k\) point index satisfying \(\mathbf{k}=\mathbf{Q}_i\mathbf{k}_j+\mathbf{G}\)./Hamiltonian/dims
: Descriptor array of length 8 containing \([0,0,0,M,N_\alpha,N_\beta,0,0]\)./Hamiltonian/H1_kp[n]
Contains the \([m_{\mathbf{k}_n},m_{\mathbf{k}_n}]\) dimensional onebody Hamiltonian matrix elements \(h_{(\mathbf{k}_{n}p)(\mathbf{k}_{n}q)}\)./Hamiltonian/ComplexIntegrals
Length 1 array that specifies if integrals are complex valued. 1 for complex integrals, 0 for real integrals./Hamiltonian/KPoints
\([N_k,3]\) Dimensional array containing \(k\)points used to sample Brillouin zone./Hamiltonian/dims
Descriptor array of length 8 containing \([0,0,N_k,M,N_\alpha,N_\beta,0,N_\mathrm{nchol}]\). Note that \(M\) is the total number of basis functions, i.e. \(M=\sum_\mathbf{k} m_\mathbf{k}\), and likewise for the number of electrons./Hamiltonian/Energies
Array containing \([E_{II}, E_{\mathrm{core}}]\). \(E_{II}\) should contain ionion repulsion energy and any additional constant terms which have to be added to the total energy (such as the electronelectron interaction Madelung contribution of \(\frac{1}{2} N \xi )\). \(E_{\mathrm{core}}\) is deprecated and not used.
Complex integrals should be written as an array with an additional dimension, e.g., a 1D array should be written as a 2D array with array_hdf5[:,0]=real(1d_array)
and array_hdf5[:,1]=imag(1d_array)
. The functions afqmctools.utils.misc.from_qmcpack_complex
and afqmctools.utils.misc.to_qmcpack_complex
can be used to transform qmcpack format to complex valued numpy arrays of the appropriate shape and vice versa.
Finally, if using external tools to generate this file format, we provide a sanity checker script in utils/afqmctools/bin/test_afqmc_input.py
which will raise errors if the format does not conform to what is being used internally.
Wavefunction File formats
AFQMC allows for two types of multideterminant trial wavefunctions: nonorthogonal multi Slater determinants (NOMSD) or SHCI/CASSCF style particlehole multi Slater determinants (PHMSD).
The file formats are described below
NOMSD
h5dump n wfn.h5
HDF5 "wfn.h5" {
FILE_CONTENTS {
group /
group /Wavefunction
group /Wavefunction/NOMSD
dataset /Wavefunction/NOMSD/Psi0_alpha
dataset /Wavefunction/NOMSD/Psi0_beta
group /Wavefunction/NOMSD/PsiT_0
dataset /Wavefunction/NOMSD/PsiT_0/data_
dataset /Wavefunction/NOMSD/PsiT_0/dims
dataset /Wavefunction/NOMSD/PsiT_0/jdata_
dataset /Wavefunction/NOMSD/PsiT_0/pointers_begin_
dataset /Wavefunction/NOMSD/PsiT_0/pointers_end_
group /Wavefunction/NOMSD/PsiT_1
dataset /Wavefunction/NOMSD/PsiT_1/data_
dataset /Wavefunction/NOMSD/PsiT_1/dims
dataset /Wavefunction/NOMSD/PsiT_1/jdata_
dataset /Wavefunction/NOMSD/PsiT_1/pointers_begin_
dataset /Wavefunction/NOMSD/PsiT_1/pointers_end_
dataset /Wavefunction/NOMSD/ci_coeffs
dataset /Wavefunction/NOMSD/dims
}
}
Note that the \(\alpha\) components of the trial wavefunction are stored under
PsiT_{2n}
and the \(\beta\) components are stored under PsiT_{2n+1}
.
/Wavefunction/NOMSD/Psi0_alpha
\([M,N_\alpha]\) dimensional array \(\alpha\) component of initial walker wavefunction./Wavefunction/NOMSD/Psi0_beta
\([M,N_\beta]\) dimensional array for \(\beta\) initial walker wavefunction./Wavefunction/NOMSD/PsiT_{2n}/data_
Array of length \(nnz\) containing nonzero elements of \(n\)th \(\alpha\) component of trial wavefunction walker wavefunction. Note the conjugate transpose of the Slater matrix is stored./Wavefunction/NOMSD/PsiT_{2n}/dims
Array of length 3 containing \([M,N_{\alpha},nnz]\) where \(nnz\) is the number of nonzero elements of this Slater matrix/Wavefunction/NOMSD/PsiT_{2n}/jdata_
CSR indices array./Wavefunction/NOMSD/PsiT_{2n}/pointers_begin_
CSR format begin index pointer array./Wavefunction/NOMSD/PsiT_{2n}/pointers_end_
CSR format end index pointer array./Wavefunction/NOMSD/ci_coeffs
\(N_D\) length array of ci coefficients. Stored as complex numbers./Wavefunction/NOMSD/dims
Integer array of length 5 containing \([M,N_\alpha,N_\beta,\) walker_type \(,N_D]\)
PHMSD
h5dump n wfn.h5
HDF5 "wfn.h5" {
FILE_CONTENTS {
group /
group /Wavefunction
group /Wavefunction/PHMSD
dataset /Wavefunction/PHMSD/Psi0_alpha
dataset /Wavefunction/PHMSD/Psi0_beta
dataset /Wavefunction/PHMSD/ci_coeffs
dataset /Wavefunction/PHMSD/dims
dataset /Wavefunction/PHMSD/occs
dataset /Wavefunction/PHMSD/type
}
}
/Wavefunction/NOMSD/Psi0_alpha
\([M,N_\alpha]\) dimensional array \(\alpha\) component of initial walker wavefunction./Wavefunction/NOMSD/Psi0_beta
\([M,N_\beta]\) dimensional array for \(\beta\) initial walker wavefunction./Wavefunction/PHMSD/ci_coeffs
\(N_D\) length array of ci coefficients. Stored as complex numbers./Wavefunction/PHMSD/dims
Integer array of length 5 containing \([M,N_\alpha,N_\beta,\) walker_type \(,N_D]\)/Wavefunction/PHMSD/occs
Integer array of length \((N_\alpha+N_\beta)*N_D\) describing the determinant occupancies. For example if \((N_\alpha=N_\beta=2)\) and \(N_D=2\), \(M=4\), and if \(\Psi_T\rangle = 0,1\rangle0,1\rangle + 0,1\rangle0,2\rangle>\) then occs = \([0, 1, 4, 5, 0, 1, 4, 6]\). Note that \(\beta\) occupancies are displacd by \(M\)./Wavefunction/PHMSD/type
integer 0/1. 1 implies trial wavefunction is written in different basis than the underlying basis used for the integrals. If so a matrix of orbital coefficients is required to be written in the NOMSD format. If 0 then assume wavefunction is in same basis as integrals.
Current Feature Implementation Status
The current status of features available in QMCPACK is as follows:
Hamiltonian 
SD 
NOMSD 
PHMSD 
Real Build 
Complex Build 

Sparse 
Yes 
Yes 
Yes 
Yes 
Yes 
Dense 
Yes 
Yes 
No 
Yes 
No 
kpoint 
Yes 
No 
No 
No 
Yes 
THC 
Yes 
No 
No 
Yes 
Yes 
Hamiltonian 
SD 
NOMSD 
PHMSD 
Real Build 
Complex Build 

Sparse 
No 
No 
No 
No 
No 
Dense 
Yes 
No 
No 
Yes 
No 
kpoint 
Yes 
No 
No 
No 
Yes 
THC 
Yes 
No 
No 
Yes 
Yes 
Advice/Useful Information
AFQMC calculations are computationally expensive and require some care to obtain reasonable performance. The following is a growing list of useful advice for new users, followed by a sample input for a large calculation.
Generate Choleskydecomposed integrals with external codes instead of the 2electron integrals directly. The generation of the Cholesky factorization is faster and consumes less memory.
Use the hybrid algorithm for walker propagation. Set steps/substeps to adequate values to reduce the number of energy evaluations. This is essential when using large multideterminant expansions.
Adjust cutoffs in the wavefunction and propagator bloxks until desired accuracy is reached. The cost of the calculation will depend on these cutoffs.
Adjust ncores/nWalkers to obtain better efficiency. Larger nWalkers will lead to more efficient linear algebra operations but will increase the time per step. Larger ncores will reduce the time per step but will reduce efficiency because of inefficiencies in the parallel implementation. For large calculations, values between 6–12 for both quantities should be reasonable, depending on architecture.
...
<Hamiltonian name="ham0" type="SparseGeneral" info="info0">
<parameter name="filename">fcidump.h5</parameter>
<parameter name="cutoff_1bar">1e6</parameter>
<parameter name="cutoff_decomposition">1e5</parameter>
</Hamiltonian>
<Wavefunction name="wfn0" type="MSD" info="info0">
<parameter name="filetype">ascii</parameter>
<parameter name="filename">wfn.dat</parameter>
</Wavefunction>
<WalkerSet name="wset0">
<parameter name="walker_type">closed</parameter>
</WalkerSet>
<Propagator name="prop0" info="info0">
<parameter name="hybrid">yes</parameter>
</Propagator>
<execute wset="wset0" ham="ham0" wfn="wfn0" prop="prop0" info="info0">
<parameter name="ncores">8</parameter>
<parameter name="timestep">0.01</parameter>
<parameter name="blocks">10000</parameter>
<parameter name="steps">10</parameter>
<parameter name="substeps">5</parameter>
<parameter name="nWalkers">8</parameter>
<parameter name="ortho">5</parameter>
</execute>
afqmc
method
parameters in AFQMCInfo
Name 
Datatype 
Values 
Default 
Description 


integer 
\(\geq 0\) 
no 
Number of molecular orbitals 

integer 
\(\geq 0\) 
no 
Number of active electrons of spinup 

integer 
\(\geq 0\) 
no 
Number of active electrons of spindown 
parameters in Hamiltonian
Name 
Datatype 
Values 
Default 
Description 


argument 
Name of 


string 
no 
Name of file with the hamiltonian 


string 
hdf5 
yes 
Native HDF5based format of QMCPACK 
parameters in Wavefunction
Name 
Datatype 
Values 
Default 
Description 


argument 
name of 


argument 
MSD, PHMSD 
no 
Linear combination of (assumed nonorthogonal) Slater determinants 

string 
ascii, hdf5 
no 
CItype multideterminant wave function 
parameters in WalkerSet
Name 
Datatype 
Values 
Default 
Description 


string 
collinear 
yes 
Request a collinear walker set. 
closed 
no 
Request a closed shell (doublyoccupied) walker set. 
parameters in Propagator
Name 
Datatype 
Values 
Default 
Description 


argument 
afqmc 
afqmc 
Type of propagator 

argument 
Name of 


string 
yes 
Use hybrid propagation algorithm. 

no 
Use local energy based propagation algorithm. 
parameters in execute
Name 
Datatype 
Values 
Default 
Description 


argument 


argument 


argument 


argument 


argument 
Name of 


integer 
\(\geq 0\) 
5 
Initial number of walkers per task group 

real 
\(> 0\) 
0.01 
Time step in 1/a.u. 

integer 
\(\geq 0\) 
100 
Number of blocks 

integer 
\(> 0\) 
1 
Number of steps within a block 

integer 
\(> 0\) 
1 
Number of substeps within a step 

integer 
\(> 0\) 
1 
Number of steps between walker orthogonalization. 
AFQMCTOOLS
The afqmctools
library found in qmcpack/utils/afqmctools
provides a number of
tools to interface electronic structure codes with AFQMC in QMCPACK. Currently PYSCF is
the best supported package and is capable of generating both molecular and solid state
input for AFQMC.
In what follows we will document the most useful routines from a user’s perspective.
afqmctools has to be in your PYTHONPATH.
pyscf_to_afqmc.py
This is the main script to convert PYSCF output into QMCPACK input. The command line options are as follows:
> pyscf_to_afqmc.py h
usage: pyscf_to_afqmc.py [h] [i CHK_FILE] [o HAMIL_FILE] [w WFN_FILE]
[q QMC_INPUT] [t THRESH] [k] [densityfit] [a]
[c CAS] [d] [n NDET_MAX] [r] [p]
[low LOW_THRESH] [high HIGH_THRESH] [dense]
[v]
optional arguments:
h, help show this help message and exit
i CHK_FILE, input CHK_FILE
Input pyscf .chk file.
o HAMIL_FILE, output HAMIL_FILE
Output file name for QMCPACK hamiltonian.
w WFN_FILE, wavefunction WFN_FILE
Output file name for QMCPACK wavefunction. By default
will write to hamil_file.
q QMC_INPUT, qmcpackinput QMC_INPUT
Generate skeleton QMCPACK input xml file.
t THRESH, choleskythreshold THRESH
Cholesky convergence threshold.
k, kpoint Generate explicit kpoint dependent integrals.
densityfit Use density fitting integrals stored in input pyscf
chkpoint file.
a, ao, orthoao Transform to ortho AO basis. Default assumes we work
in MO basis
c CAS, cas CAS Specify a CAS in the form of N,M.
d, disableham Disable hamiltonian generation.
n NDET_MAX, numdets NDET_MAX
Set upper limit on number of determinants to generate.
r, realham Write integrals as real numbers.
p, phdf Use parallel hdf5.
low LOW_THRESH Lower threshold for noninteger occupanciesto include
in multideterminant exansion.
high HIGH_THRESH Upper threshold for noninteger occupanciesto include
in multideterminant exansion.
dense Write dense Hamiltonian.
v, verbose Verbose output.
examples on how to generate AFQMC input from PYSCF simulations are available in AFQMC Tutorials
afqmc_to_fcidump.py
This script is useful for converting AFQMC hamiltonians to the FCIDUMP format.
> afqmc_to_fcidump.py
usage: afqmc_to_fcidump.py [h] [i INPUT_FILE] [o OUTPUT_FILE] [s SYMM]
[t TOL] [c] [complexparen] [v]
optional arguments:
h, help show this help message and exit
i INPUT_FILE, input INPUT_FILE
Input AFQMC hamiltonian file.
o OUTPUT_FILE, output OUTPUT_FILE
Output file for FCIDUMP.
s SYMM, symmetry SYMM
Symmetry of integral file (1,4,8).
t TOL, tol TOL Cutoff for integrals.
c, complex Whether to write integrals as complex numbers.
complexparen Whether to write FORTRAN format complex numbers.
v, verbose Verbose output.
fcidump_to_afqmc.py
This script is useful for converting Hamiltonians in the FCIDUMP format to the AFQMC file format.
> fcidump_to_afqmc.py h
usage: fcidump_to_afqmc.py [h] [i INPUT_FILE] [o OUTPUT_FILE]
[writecomplex] [t THRESH] [s SYMM] [v]
optional arguments:
h, help show this help message and exit
i INPUT_FILE, input INPUT_FILE
Input FCIDUMP file.
o OUTPUT_FILE, output OUTPUT_FILE
Output file name for PAUXY data.
writecomplex Output integrals in complex format.
t THRESH, choleskythreshold THRESH
Cholesky convergence threshold.
s SYMM, symmetry SYMM
Symmetry of integral file (1,4,8).
v, verbose Verbose output.
Writing a Hamiltonian
write_qmcpack_sparse
and write_qmcpack_dense
can be used to write either sparse or
dense qmcpack Hamiltonians.
import numpy
from afqmctools.hamiltonian.io import write_qmcpack_sparse, write_qmcpack_dense
nmo = 50
nchol = 37
nelec = (3,3)
enuc = 108.3
# hcore and eri should obey the proper symmetry in real applications
# h_ij
hcore = numpy.random.random((nmo,nmo))
# L_{(ik),n}
chol = numpy.random.random((nmo*nmo, nchol))
write_qmcpack_dense(hcore, chol, nelec, nmo, enuc,
real_chol=True,
filename='hamil_dense.h5')
write_qmcpack_sparse(hcore, chol, nelec, nmo, enuc,
real_chol=True,
filename='hamil_sparse.h5')
Note the real_chol
parameter controls whether the integrals are written as real or
complex numbers. Complex numbers should be used if DENABLE_QMC_COMPLEX=1
, while the
dense Hamiltonian is only available for real builds.
Writing a wavefunction
write_qmcpack_wfn
can be used to write either NOMSD or PHMSD wavefunctions:
import numpy
from afqmctools.wavefunction.mol import write_qmcpack_wfn
# NOMSD
ndet = 100
nmo = 50
nelec = (3, 7)
wfn = numpy.array(numpy.random.random((ndet, nmo, sum(nelec))), dtype=numpy.complex128)
coeffs = numpy.array(numpy.random.random((ndet)), dtype=numpy.complex128)
uhf = True
write_qmcpack_wfn('wfn.h5', (coeffs, wfn), uhf, nelec, nmo)
By default the first term in the expansion will be used as the initial walker
wavefunction. To use another wavefunction we can pass a value to the init
parameter:
init = numpy.array(numpy.random.random((nmo,sum(nelec)), dtype=numpy.complex128)
write_qmcpack_wfn('wfn.h5', (coeffs, wfn), uhf, nelec, nmo, init=[init,init])
Particlehole wavefunction (PHMSD) from SHCI or CASSCF calculations are also written using the same function:
import numpy
from afqmctools.wavefunction.mol import write_qmcpack_wfn
# PHMSD
ndet = 2
nmo = 4
nelec = (2,2)
uhf = True
# psi_T> = 1/sqrt(2)(0,1>0,1> + 0,1>0,2>)
coeffs = numpy.array([0.707,0.707], dtype=numpy.complex128)
occa = numpy.array([(0,1), (0,1)])
occb = numpy.array([(0,1), (0,2)])
write_qmcpack_wfn('wfn.h5', (coeffs, occa, occb), uhf, nelec, nmo)
Analyzing Estimators
The afqmctools.analysis.average
module can be used to perform simple error analysis
for estimators computed with AFQMC.
Warning
Autocorrelation is not accounted for. Use with caution.
 average_one_rdm
Returns P[s,i,j] = \(\langle c_{is}^{\dagger} c_{js}\rangle\) as a (nspin, M, M) dimensional array.
 average_two_rdm
Gamma[s1s2,i,k,j,l] = \(\langle c_{i}^{\dagger} c_{j}^{\dagger} c_{l} c_{k} \rangle\). For closed shell systems, returns [(a,a,a,a),(a,a,b,b)]. For collinear systems, returns [(a,a,a,a),(a,a,b,b),(b,b,b,b)].
 average_diag_two_rdm
Returns \(\langle c_{is}^+ c_{jt}^+ c_{jt} c_{is}\rangle\) as a (2M,2M) dimensional array.
 average_on_top_pdm
Returns \(n_2(\mathbf{r},\mathbf{r})\) for a given real space grid.
 average_realspace_correlations
Returns \(\langle C(\mathbf{r}_1)C(\mathbf{r}_2) \rangle\) and \(\langle S(\mathbf{r}_1)S(\mathbf{r}_2) \rangle\) for a given set of points in real space. \(\hat{C} = (\hat{n}_\uparrow+ \hat{n}_\downarrow)\), \(\hat{S}=(\hat{n}_\uparrow\hat{n}_\downarrow)\)
 average_atom_correlations
Returns \(\langle C(I) \rangle\), \(\langle S(I) \rangle\), \(\langle C(I) C(J) \rangle\), \(\langle S(I) S(J) \rangle\) for a given set of atomic sites \(I,J\). \(\hat{C} = (\hat{n}_\uparrow+ \hat{n}_\downarrow)\), \(\hat{S}=(\hat{n}_\uparrow\hat{n}_\downarrow)\)
 average_gen_fock
Returns generalized Fock matrix \(F_{\pm}\). The parameter
fock_type
is used to specify \(F_{+}\) (fock_type='plus'
) or \(F_{}\) (fock_type='minus'
) get_noons
Get natural orbital occupation numbers from onerdm.
As an example the following will extract the back propagated one rdm for the maximum propagation time, and skip 10 blocks as the equilibration phase.
from afqmctools.analysis.average import average_one_rdm
P, Perr = average_one_rdm('qmc.s000.stat.h5', estimator='back_propagated', eqlb=10)
 ADVFerre+09
Francesco Aquilante, Luca De Vico, Nicolas Ferré, Giovanni Ghigo, Peråke Malmqvist, Pavel Neogrády, Thomas Bondo Pedersen, Michal Pitoňák, Markus Reiher, Björn O. Roos, Luis SerranoAndrés, Miroslav Urban, Valera Veryazov, and Roland Lindh. MOLCAS 7: The Next Generation. J. Comput. Chem., 31(1):224, 2009. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21318.
 BL77
Nelson H. F. Beebe and Jan Linderberg. Simplifications in the generation and transformation of twoelectron integrals in molecular calculations. Int. J. Quantum Chem., 12(4):683, 1977. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.560120408.
 HPMartinez12
Edward G. Hohenstein, Robert M. Parrish, and Todd J. Martínez. Tensor hypercontraction density fitting. I. Quartic scaling second and thirdorder MøllerPlesset perturbation theory. The Journal of Chemical Physics, 137(4):044103, 2012. URL: https://doi.org/10.1063/1.4732310, doi:10.1063/1.4732310.
 HPSMartinez12
Edward G. Hohenstein, Robert M. Parrish, C. David Sherrill, and Todd J. Martínez. Communication: Tensor hypercontraction. III. Leastsquares tensor hypercontraction for the determination of correlated wavefunctions. The Journal of Chemical Physics, 137(22):221101, 2012. URL: https://doi.org/10.1063/1.4768241, doi:10.1063/1.4768241.
 KdMerasP03
Henrik Koch, Alfredo Sánchez de Merás, and Thomas Bondo Pedersen. Reduced scaling in electronic structure calculations using Cholesky decompositions. The Journal of Chemical Physics, 118(21):9481, 2003. URL: https://doi.org/10.1063/1.1578621.
 MZM20
Fionn D Malone, Shuai Zhang, and Miguel A Morales. Accelerating auxiliaryfield quantum monte carlo simulations of solids with graphical processing unit. arXiv preprint arXiv:2003.09468, 2020.
 MZM19
Fionn D. Malone, Shuai Zhang, and Miguel A. Morales. Overcoming the Memory Bottleneck in Auxiliary Field Quantum Monte Carlo Simulations with Interpolative Separable Density Fitting. J. Chem. Theory. Comput., 15(1):256, 2019. URL: https://doi.org/10.1021/acs.jctc.8b00944, doi:10.1021/acs.jctc.8b00944.
 MZC19
Mario Motta, Shiwei Zhang, and Garnet KinLic Chan. Hamiltonian symmetries in auxiliaryfield quantum Monte Carlo calculations for electronic structure. Physical Review B, 100:045127, July 2019. URL: https://link.aps.org/doi/10.1103/PhysRevB.100.045127, doi:10.1103/PhysRevB.100.045127.
 PHMartinezS12
Robert M. Parrish, Edward G. Hohenstein, Todd J. Martínez, and C. David Sherrill. Tensor hypercontraction. II. Leastsquares renormalization. The Journal of Chemical Physics, 137(22):224106, 2012. URL: https://doi.org/10.1063/1.4768233, doi:10.1063/1.4768233.
 PKVZ11
Wirawan Purwanto, Henry Krakauer, Yudistira Virgus, and Shiwei Zhang. Assessing weak hydrogen binding on Ca+ centers: An accurate manybody study with large basis sets. The Journal of Chemical Physics, 135(16):164105, 2011. URL: https://doi.org/10.1063/1.3654002, doi:10.1063/1.3654002.
 PZ04
Wirawan Purwanto and Shiwei Zhang. Quantum monte carlo method for the ground state of manyboson systems. Phys. Rev. E, 70:056702, November 2004. doi:10.1103/PhysRevE.70.056702.
 PZK13
Wirawan Purwanto, Shiwei Zhang, and Henry Krakauer. FrozenOrbital and Downfolding Calculations with AuxiliaryField Quantum Monte Carlo. Journal of Chemical Theory and Computation, 9(11):4825–4833, 2013. URL: https://doi.org/10.1021/ct4006486.
 Zha13
Shiwei Zhang. Auxiliaryfield quantum monte carlo for correlated electron systems. Modeling and Simulation, 3:, 2013. URL: http://hdl.handle.net/2128/5389, doi:.
 ZK03(1,2)
Shiwei Zhang and Henry Krakauer. Quantum monte carlo method using phasefree random walks with slater determinants. Phys. Rev. Lett., 90:136401, April 2003. doi:10.1103/PhysRevLett.90.136401.