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

Merge pull request #1755 from prckent/improve_chapters 

Further manual editing
This commit is contained in:
Paul R. C. Kent 2019-07-22 13:10:06 -04:00 committed by GitHub
parent 2dacbce6e8 2aa6f04ced
commit 04e4ee60b9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 605 additions and 618 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
manual/afqmc.tex
View File

@ -1,14 +1,14 @@
\chapter{Auxiliary-Field Quantum Monte Carlo}
\label{chap:afqmc}
The Auxiliary-Field Quantum Monte Carlo (AFQMC) method is an orbital-space formulation of the imaginary-time propagation algorithm. We refer the reader to one of the review articles on the method \cite{AFQMC_review,PhysRevLett.90.136401,PhysRevE.70.056702}, for a detailed description of the algorithm. It uses the Hubbard-Stratonovich transformation to express the imaginary-time propagator, which is inherently a 2-body operator, as an integral over 1-body propagators which can be efficiently applied to an arbitrary Slater determinant. This transformation allows us to represent the interacting many-body system as an average over a non-interacting system (e.g. Slater determinants) in a time-dependent fluctuating external field (the Auxiliary fields). The walkers in this case represent non-orthogonal Slater determinants, whose time average represent the desired quantum state. QMCPACK currently implements the phaseless AFQMC algorithm of Zhang and Krakauer \cite{PhysRevLett.90.136401}, where a trial wave-function 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 fixed-node approximation in real-space DMC, but applied in the Hilbert space where the AFQMC random walk occurs.
The AFQMC method is an orbital-space formulation of the imaginary-time propagation algorithm. We refer the reader to one of the review articles on the method \cite{AFQMC_review,PhysRevLett.90.136401,PhysRevE.70.056702} for a detailed description of the algorithm. It uses the Hubbard-Stratonovich transformation to express the imaginary-time propagator, which is inherently a 2-body operator, as an integral over 1-body propagators, which can be efficiently applied to an arbitrary Slater determinant. This transformation allows us to represent the interacting many-body system as an average over a noninteracting system (e.g., Slater determinants) in a time-dependent 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 \cite{PhysRevLett.90.136401}, 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 fixed-node approximation in real-space DMC but applied in the Hilbert space where the AFQMC random walk occurs.
\section{Theoretical Background}
... Coming Soon ...
\section{Input}
The input for an AFQMC calculation is fundamentally different to the input for other real-space 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 are found below. The input file has six basic xml-blocks: \texttt{AFQMCInfo}, \texttt{Hamiltonian}, \texttt{Wavefunction}, \texttt{WalkerSet}, \texttt{Propagator}, and \texttt{execute}. The first five define input structures required for various types of calculations. The \texttt{execute} block represents actual calculations and takes as input the other blocks.
Non-execution blocks are parsed first, followed by a second pass where execution blocks are parsed (and executed) in order. Listing 13.1 shows an example of a minimal input file for an AFQMC calculation. Table~\ref{table:afqmc_basic} 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.
The input for an AFQMC calculation is fundamentally different to the input for other real-space 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 xml-blocks: \texttt{AFQMCInfo}, \texttt{Hamiltonian}, \texttt{Wavefunction}, \texttt{WalkerSet}, \texttt{Propagator}, and \texttt{execute}. The first five define input structures required for various types of calculations. The \texttt{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 15.1 shows an example of a minimal input file for an AFQMC calculation. Table \ref{table:afqmc_basic} 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.
\begin{lstlisting}[style=QMCPXML,caption=Sample input file for AFQMC.]
<?xml version="1.0"?>
@ -50,22 +50,23 @@ Non-execution blocks are parsed first, followed by a second pass where execution
%The following table lists some of the practical parameters
\begin{table}[h]
\begin{center}
\caption{Input options for AFQMC in QMCPACK.\label{table:afqmc_basic}}
\begin{tabularx}{\textwidth}{l l l l l X }
\hline
\multicolumn{6}{l}{\texttt{afqmc} method} \\
\hline
\multicolumn{6}{l}{parameters in \texttt{AFQMCInfo}} \\
& \bfseries name & \bfseries datatype & \bfseries values & \bfseries default & \bfseries description \\
& \texttt{NMO } & integer & $\ge 0$ & no & number of molecular orbitals \\
& \texttt{NAEA } & integer & $\ge 0$ & no & number of active electrons of spin up \\
& \texttt{NAEB } & integer & $\ge 0$ & no & number of active electrons of spin down \\
& \texttt{NMO } & integer & $\ge 0$ & no & Number of molecular orbitals \\
& \texttt{NAEA } & integer & $\ge 0$ & no & Number of active electrons of spin-up \\
& \texttt{NAEB } & integer & $\ge 0$ & no & Number of active electrons of spin-down \\
\multicolumn{6}{l}{parameters in \texttt{Hamiltonian}} \\
& \texttt{info } & argument & & & name of \texttt{AFQMCInfo} block \\\\
& \texttt{filename } & string & & no & name of file with the hamiltonian \\
& \texttt{filetype } & string & hdf5 & yes & native HDF5-based format of QMCPACK \\
& \texttt{info } & argument & & & Name of \texttt{AFQMCInfo} block \\\\
& \texttt{filename } & string & & no & Name of file with the hamiltonian \\
& \texttt{filetype } & string & hdf5 & yes & Native HDF5-based format of QMCPACK \\
\multicolumn{6}{l}{parameters in \texttt{Wavefunction}}\\
& \texttt{info } & argument & & & name of \texttt{AFQMCInfo} block \\
& \texttt{type } & argument & MSD & no & linear combination of (assumed non-orthogonal) Slater determinants \\
& \texttt{info } & argument & & & Name of \texttt{AFQMCInfo} block \\
& \texttt{type } & argument & MSD & no & Linear combination of (assumed non-orthogonal) Slater determinants \\
& \texttt{ } & & PHMSD & & CI-type multi-determinant wave function \\
& \texttt{filetype } & string & ascii & no & ASCII data file type \\
& \texttt{ } & & hdf5 & & HDF5 data file type \\
@ -73,8 +74,8 @@ Non-execution blocks are parsed first, followed by a second pass where execution
& \texttt{walker$\_$type } & string & collinear & yes & Request a collinear walker set. \\
& \texttt{ } & & closed & no & Request a closed shell (doubly-occupied) walker set. \\
\multicolumn{6}{l}{parameters in \texttt{Propagator}} \\
& \texttt{type } & argument & afqmc & afqmc & type of propagator \\
& \texttt{info } & argument & & & name of \texttt{AFQMCInfo} block \\
& \texttt{type } & argument & afqmc & afqmc & Type of propagator \\
& \texttt{info } & argument & & & Name of \texttt{AFQMCInfo} block \\
& \texttt{hybrid } & string & yes & yes & Use hybrid propagation algorithm. \\
& \texttt{ } & & no & & Use local energy based propagation algorithm. \\
\multicolumn{6}{l}{parameters in \texttt{execute}} \\
@ -82,51 +83,49 @@ Non-execution blocks are parsed first, followed by a second pass where execution
& \texttt{ham } & argument & & & \\
& \texttt{wfn } & argument & & & \\
& \texttt{prop } & argument & & & \\
& \texttt{info } & argument & & & name of \texttt{AFQMCInfo} block \\
& \texttt{nWalkers } & integer & $\ge 0$ & 5 & initial number of walkers per task group \\
& \texttt{timestep } & real & $> 0$ & 0.01 & time step in 1/a.u. \\
& \texttt{blocks } & integer & $\ge 0$ & 100 & number of blocks \\
& \texttt{step } & integer & $> 0$ & 1 & number of steps within a block \\
& \texttt{substep } & integer & $> 0$ & 1 & number of substeps within a step \\
& \texttt{ortho } & integer & $> 0$ & 1 & number of steps between walker orthogonalization. \\
& \texttt{info } & argument & & & Name of \texttt{AFQMCInfo} block \\
& \texttt{nWalkers } & integer & $\ge 0$ & 5 & Initial number of walkers per task group \\
& \texttt{timestep } & real & $> 0$ & 0.01 & Time step in 1/a.u. \\
& \texttt{blocks } & integer & $\ge 0$ & 100 & Number of blocks \\
& \texttt{step } & integer & $> 0$ & 1 & Number of steps within a block \\
& \texttt{substep } & integer & $> 0$ & 1 & Number of substeps within a step \\
& \texttt{ortho } & integer & $> 0$ & 1 & Number of steps between walker orthogonalization. \\
\hline
\end{tabularx}
\end{center}
\caption{Input options for AFQMC in QMCPACK}
\label{table:afqmc_basic}
\end{table}
Below is a list of 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 can change in the future.\\
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.\\
\texttt{AFQMCInfo}: input block that defines basic information about the calculation. It is passed to all other input blocks to propagate the basic information:
\texttt{AFQMCInfo}: Input block that defines basic information about the calculation. It is passed to all other input blocks to propagate the basic information:
\texttt{<AFQMCInfo name="info0">}
\begin{itemize}
\item \textbf{NMO}. Number of molecular orbitals, i.e., number of states in the single particle basis.
\item \textbf{NAEA}. Number of Active Electrons-Alpha, i.e., number of spin-up electrons.
\item \textbf{NAEB}. Number of Active Electrons-Beta, i.e., number of spin-down electrons.
\item \textbf{NAEA}. Number of active electrons-alpha, i.e., number of spin-up electrons.
\item \textbf{NAEB}. Number of active electrons-beta, i.e., number of spin-down electrons.
\end{itemize}
\texttt{Hamiltonian}: controls the object that reads, stores and manages the hamiltonian.
\texttt{Hamiltonian}: Controls the object that reads, stores, and manages the \texttt{hamiltonian}.
\texttt{<Hamiltonian name="ham0" type="SparseGeneral" info="info0">}
\begin{itemize}
\item \textbf{filename}. Name of file with the \texttt{Hamiltonian}. This is a required parameter.
\item \textbf{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 only applied to the 2-electron integrals). Default: 1e-8
\item \textbf{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 2-electron integrals). Default: 1e-8
\item \textbf{cutoff\_decomposition}. Cutoff used to stop the iterative cycle in the generation of the Cholesky decomposition of the 2-electron 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: 1e-6
\item \textbf{nblocks}. This parameter controls the distribution of the 2-electron 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
\item \textbf{printEig}. If ``yes'', prints additional information during the Cholesky decomposition.
Default: no
\item \textbf{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 2-electron integrals are not positive definite due to round-off errors in their generation.
\item \textbf{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 2-electron integrals are not positive definite because of round-off errors in their generation.
Default: no \\
\end{itemize}
\texttt{Wavefunction}: controls the object that manages the trial wave-functions. This block expects a list of xml-blocks defining actual trial-wave functions for various roles.
\texttt{Wavefunction}: controls the object that manages the trial wavefunctions. This block expects a list of xml-blocks defining actual trial wavefunctions for various roles.
\texttt{<Wavefunction name="wfn0" type="MSD/PHMSD" info="info0">}
\begin{itemize}
\item \textbf{filename}. Name of file with wave-function information.
\item \textbf{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.
\item \textbf{filename}. Name of file with wavefunction information.
\item \textbf{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: 1e-6
\item \textbf{nnodes}. Defines the parallelization of the local energy evaluation and the distribution of the \texttt{Hamiltonian} matrix (not to be confused with the list of 2-electron integrals managed by \texttt{Hamiltonian}. These are not the same.) If nnodes $>$ 1, the nodes in the simulation are split into groups of nnodes, each group works collectively in the evaluation of the local energy of their walkers. This helps distribute the effort involved in the evaluation of the local energy among the nodes in the group, but also distributes the memory associated with the wave-function among the nodes in the group.
\item \textbf{nnodes}. Defines the parallelization of the local energy evaluation and the distribution of the \texttt{Hamiltonian} matrix (not to be confused with the list of 2-electron integrals managed by \texttt{Hamiltonian}. These are not the same.) If nnodes $>$ 1, the nodes in the simulation are split into groups of nnodes, each group works collectively in the evaluation of the local energy of their walkers. This helps distribute the effort involved in the evaluation of the local energy among the nodes in the group, but also distributes the memory associated with the wavefunction among the nodes in the group.
Default: No distribution
\item \textbf{ndet}. Number of determinants to read from file.
Default: Read all determinants.
@ -154,13 +153,13 @@ Below is a list of all input sections for AFQMC calculations, along with a detai
\begin{itemize}
\item \textbf{cutoff}. Cutoff applied to Cholesky vectors. Elements of the Cholesky vectors below this value are set to zero. Only meaningful with sparse hamiltonians.
Default: 1e-6
\item \textbf{substractMF}. If ``yes'', apply mean-field subtraction based on the ImpSamp trial wave-function. Must set to ``no'' to turn it off.
\item \textbf{substractMF}. If ``yes'', apply mean-field subtraction based on the ImpSamp trial wavefunction. Must set to ``no'' to turn it off.
Default: yes
\item \textbf{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 wave-function).
\item \textbf{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
\item \textbf{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
\item \textbf{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 wave-function is never evaluated. To obtain energy estimates in this case, you must define an Estimator xml-block with the \texttt{Wavefunction} block. The local energy of this trial wave-function is evaluated and printed. It is possible to use a previously defined trial wave-function in the Estimator block, just set its ``name'' argument to the name of a previously defined wave-function. In this case, the same object is used for both roles.
\item \textbf{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 xml-block with the \texttt{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: no
\item \textbf{nnodes}. 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)
@ -173,17 +172,17 @@ Below is a list of all input sections for AFQMC calculations, along with a detai
\texttt{execute}: Defines an execution region.
\texttt{<execute wset="wset0" ham="ham0" wfn="wfn0" prop="prop0" info="info0">}
\begin{itemize}
\item \textbf{nWalkers}. initial number of walkers per core group (see ncores). This sets the number of walkers for a given gorup 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 $. \\
\item \textbf{nWalkers}. Initial number of walkers per core group (see ncores). This sets the number of walkers for a given gorup 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
\item \textbf{timestep}. time step in 1/a.u. \\
\item \textbf{timestep}. Time step in 1/a.u. \\
Default: 0.01
\item \textbf{blocks}. number of blocks. Slow operations occur once per block, e.g. write to file, slow observables, checkpoints, etc. \\
\item \textbf{blocks}. Number of blocks. Slow operations occur once per block (e.g., write to file, slow observables, checkpoints), \\
Default: 100
\item \textbf{step}. number of steps within a block. Operations that occur at the step level include: load balance, orthogonalization, branching, etc. \\
\item \textbf{step}. Number of steps within a block. Operations that occur at the step level include load balance, orthogonalization, branching, etc. \\
Default: 1
\item \textbf{substep}. number of substeps within a step. Only walker propagation occurs in a substep. \\
\item \textbf{substep}. Number of substeps within a step. Only walker propagation occurs in a substep. \\
Default: 1
\item \textbf{ortho}. number of steps between orthogonalization.
\item \textbf{ortho}. Number of steps between orthogonalization.
Default: 1
\item \textbf{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 \texttt{Wavefunction} and \texttt{Propagator} task groups. The walker sets are shares by the ncores on a given node in the task group.
\item \textbf{checkpoint}. Number of blocks between checkpoint files are generated. If a value smaller than 1 is given, no file is generated. If \textbf{hdf\_write\_file} is not set, a default name is used. \textbf{Default: 0}
@ -218,13 +217,13 @@ The back\_propagation estimator has the following parameters:
\section{Advice/Useful Information}
AFQMC calculations are computationally expensive and require some care in order to obtain reasonable performance.
Below is a growing list of useful advice for new users followed by a sample input for a large calculation.
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.
\begin{itemize}
\item Generate Cholesky-decomposed integrals with external codes instead of the 2-electron integrals directly. The generation of the Cholesky factorization is faster and consumes less memory.
\item 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 multi-determinant expansions.
\item Adjust cutoffs in the wave-function and propagator bloxks until desired accuracy is reached. The cost of the calculation will depend on these cutoffs.
\item 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 due to inefficiencies in the parallel implementation. For large calculations, values between 6-12 for both quantities should be reasonable, depending on architecture.
\item 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.
\item Adjust cutoffs in the wavefunction and propagator bloxks until desired accuracy is reached. The cost of the calculation will depend on these cutoffs.
\item 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.
\end{itemize}
\begin{lstlisting}[style=QMCPXML,caption=Example of sections of an AFQMC input file for a large calculation.]

8
manual/examples.tex
View File

@ -7,7 +7,7 @@ The following examples should run in serial on a modern workstation in a few hou
\section{Using QMCPACK directly}
In \ishell{examples/molecules}, there are the following examples.
In \ishell{examples/molecules} are the following examples.
Each directory also contains a \ishell{README} file with more details.
\begin{tabular}{l l}
@ -24,7 +24,7 @@ Directory & Description \\
For more information about Nexus, see the User Guide in \ishell{nexus/documentation}.
For Python to find the Nexus library, the PYTHONPATH environment variable should be set to \ishell{<QMCPACK source>/nexus/library}.
For these examples to work properly, the executables for Quantum ESPRESSO and QMCPACK either
For these examples to work properly, the executables for QE and QMCPACK either
need to be on the path, or the paths in the script should be adjusted.
These examples can be found under the \ishell{nexus/examples/qmcpack} directory.
@ -43,8 +43,8 @@ Directory & Description \\
\ishell{graphene} & Graphene sheet with DMC \\
\ishell{c20} & C20 cage molecule \\
\ishell{oxygen\_dimer} & Binding curve for O$_2$ molecule \\
\ishell{H2O} & H$_2$O molecule with Quantum ESPRESSO orbitals \\
\ishell{LiH} & LiH crystal with Quantum ESPRESSO orbitals \\
\ishell{H2O} & H$_2$O molecule with QE orbitals \\
\ishell{LiH} & LiH crystal with QE orbitals \\
\end{tabular}

571
manual/lab_advanced_molecules.tex
View File

@ -1,15 +1,15 @@
\chapter{Lab 3: Advanced Molecular Calculations}
\chapter{Lab 3: Advanced molecular calculations}
\label{chap:lab_advanced_molecules}
\section{Topics covered in this Lab}
\section{Topics covered in this lab}
This lab covers molecular QMC calculations with wavefunctions of increasing sophistication. All of the trial wavefunctions are initially generated with the GAMESS code. Topics covered include:
\begin{itemize}
\item{Generating single determinant trial wavefunctions with GAMESS (HF and DFT)}
\item{Generating multi-determinant trial wavefunctions with GAMESS (CISD, CASCI, SOCI)}
\item{Generating single-determinant trial wavefunctions with GAMESS (HF and DFT)}
\item{Generating multideterminant trial wavefunctions with GAMESS (CISD, CASCI, and SOCI)}
\item{Optimizing wavefunctions (Jastrow factors and CSF coefficients) with QMC}
\item{DMC timestep and walker population convergence studies}
\item{DMC time step and walker population convergence studies}
\item{Systematic progressions of Jastrow factors in VMC}
\item{Systematic convergence of DMC energies with multi-determinant wavefunctions}
\item{Systematic convergence of DMC energies with multideterminant wavefunctions}
\item{Influence of orbitals basis choice on DMC energy}
%\item{Influence of orbitals basis choice on rate of multi-determinant DMC convergence}
\end{itemize}
@ -73,25 +73,25 @@ labs/lab3_advanced_molecules/exercises
\section{Exercise \#1: Basics}
The purpose of this exercise is to show how to generate wave-functions for QMCPACK
using GAMESS and to optimize the resulting wave-functions using VMC. This will be
followed by a study of the time-step and walker population dependence of DMC energies.
The purpose of this exercise is to show how to generate wavefunctions for QMCPACK
using GAMESS and to optimize the resulting wavefunctions using VMC. This will be
followed by a study of the time step and walker population dependence of DMC energies.
The exercise will be performed on a water molecule at the equilibrium geometry.
\subsection{Generation of a Hatree-Fock wave-function with GAMESS}
\subsection{Generation of a Hartree-Fock wavefunction with GAMESS}
From the top directory, go to ``ex1\_first-run-hartree-fock/gms''. This directory contains an input
From the top directory, go to ``\texttt{ex1\_first-run-hartree-fock/gms}.'' This directory contains an input
file for a HF calculation of a water molecule using BFD ECPs and the corresponding
cc-pVTZ basis set. The input file should be named: “h2o.hf.inp”. Study the input
file. If the student wishes, he can refer to section A for a more detailed description of the
GAMESS input syntax. There will be a better time to do this soon, so we recommend
that the student continues with the exercise at this point. After you are done, execute
GAMESS with this input and store the standard output in a file named ``h2o.hf.output''.
Finally, in the ``convert'' folder, use convert4qmc to generate the QMCPACK particleset and wavefunction files. It
is always useful to rename the files generated by convert4qmc to something meaningful,
since by default they are called sample.Gaussian-G2.xml and sample.Gaussian-G2.ptcl.xml.
In a standard computer (without cross-compilation), these tasks could be accomplished by
cc-pVTZ basis set. The input file should be named: “h2o.hf.inp.” Study the input
file. See Section~\ref{sec:lab_adv_mol_gamess}, ``Appendix A: GAMESS input" for a more detailed description of the
GAMESS input syntax. However, there will be a better time to do this soon, so we recommend
continuing with the exercise at this point. After you are done, execute
GAMESS with this input and store the standard output in a file named ``h2o.hf.output.''
Finally, in the ``convert'' folder, use \texttt{convert4qmc} to generate the QMCPACK \texttt{particleset} and \texttt{wavefunction} files. It
is always useful to rename the files generated by \texttt{convert4qmc} to something meaningful
since by default they are called \texttt{sample.Gaussian-G2.xml} and \texttt{sample.Gaussian-G2.ptcl.xml}.
In a standard computer (without cross-compilation), these tasks can be accomplished by
the following commands.
\begin{lstlisting}[style=SHELL]
cd ${TRAINING TOP}/ex1_first-run-hartree-fock/gms
@ -124,29 +124,29 @@ grep "TOTAL ENERGY =" h2o.hf.output
%submit the execution of the converter using ``qsub submit convert.csh'' (all QMCPACK execution
%scripts will be submitted with qsub). This is a good time to review section B, which
%contains a description on the use of the converter.
As the job runs on VESTA, it is a good time to review section B, which contains a description on the use of the converter.
As the job runs on VESTA, it is a good time to review Section~\ref{sec:lab_adv_mol_convert4qmc}, ``Appendix B: convert4qmc," which contains a description on the use of the converter.
\subsection{Optimize the wave-function}
When the execution of the previous steps is completed, there should be 2 new
files called h2o.wfs.xml and h2o.ptcl.xml. Now we will use VMC to optimize the
Jastrow parameters in the wave-function. From the top directory, go to
``ex1\_first-run-hartree-fock/opt''. Copy the xml files generated in the previous step
\subsection{Optimize the wavefunction}
When execution of the previous steps is completed, there should be two new
files called \texttt{h2o.wfs.xml} and \texttt{h2o.ptcl.xml}. Now we will use VMC to optimize the
Jastrow parameters in the wavefunction. From the top directory, go to
``\texttt{ex1\_first-run-hartree-fock/opt}.'' Copy the xml files generated in the previous step
to the current directory. This directory should already contain a basic QMCPACK input
file for an optimization calculation (optm.xml) %and a submission script (submit.csh).
Open optm.xml with your favorite text editor and modify the name of the files that contain the
wavefunction and particleset XML blocks. These files are included with the commands:
file for an optimization calculation (\texttt{optm.xml}) %and a submission script (submit.csh).
Open \texttt{optm.xml} with your favorite text editor and modify the name of the files that contain the
\texttt{wavefunction} and \texttt{particleset} XML blocks. These files are included with the commands:
\begin{lstlisting}[style=QMCPXML]
<include href=ptcl.xml/>
<include href=wfs.xml/>
\end{lstlisting}
(the particle set must be
defined before the wave-function). The name of the particle set and wave-function files should now be h2o.ptcl.xml
and h2o.wfs.xml, respectively. Study both files and submit when you are ready. Notice that the
location of the ECPs has been set for you, in your own calculations you have to make
defined before the wavefunction). The name of the particle set and wavefunction files should now be \texttt{h2o.ptcl.xml}
and \texttt{h2o.wfs.xml}, respectively. Study both files and submit when you are ready. Notice that the
location of the ECPs has been set for you; in your own calculations you have to make
sure you obtain the ECPs from the appropriate libraries and convert them to QMCPACK
format using ppconvert. This is a good time to study section C, which contains a review
of the main parameters in the optimization XML block, while this calculation finishes. The
format using ppconvert. While these calculations finish is a good time to study Section~\ref{sec:lab_adv_mol_opt_appendix}, ``Appendix C: Wavefunction optimization XML block," which contains a review
of the main parameters in the optimization XML block. The
previous steps can be accomplished by the following commands:
\begin{shade}
cd ${TRAINING TOP}/ex1_first-run-hartree-fock/opt
@ -156,9 +156,9 @@ cp ../convert/h2o.ptcl.xml ./
jobrun_vesta qmcpack optm.xml
\end{shade}
Use the analysis tool qmca to analyze the results of the calculation. Obtain the VMC
Use the analysis tool \texttt{qmca} to analyze the results of the calculation. Obtain the VMC
energy and variance for each step in the optimization and plot it using your favorite program.
Remember that qmca has built-in functions to plot the analyzed data.
Remember that \texttt{qmca} has built-in functions to plot the analyzed data.
\begin{shade}
qmca -q e *scalar.dat -p
\end{shade}
@ -171,25 +171,25 @@ qmca -q e *scalar.dat -p
\label{fig:lam_opt_conv}
\end{figure}
The resulting energy as a function of optimization step should look qualitatively similar to figure~\ref{fig:lam_opt_conv}.
The energy should decrease quickly as a function of the number of optimization steps. After 6-8 steps, the energy should be converged to $\sim$2-3mHa. To improve convergence,
we would need to increase the number of samples used during the optimization. You can
check this for yourself on your free time. With optimized wave-functions we are in a position
to perform VMC and DMC calculations. The modified wave-function files after each step
are written in a file named ID.sNNN.opt.xml, where ID is the identifier of the calculation
The resulting energy as a function of the optimization step should look qualitatively similar to Figure \ref{fig:lam_opt_conv}.
The energy should decrease quickly as a function of the number of optimization steps. After 6--8 steps, the energy should be converged to $\sim$2--3 mHa. To improve convergence,
we would need to increase the number of samples used during optimization (You can
check this for yourself later.). With optimized wavefunctions, we are in a position
to perform VMC and DMC calculations. The modified wavefunction files after each step
are written in a file named \texttt{ID.sNNN.opt.xml}, where ID is the identifier of the calculation
defined in the input file (this is defined in the project XML block with parameter “id”) and
NNN is a series number which increases with every executable xml block in the input file.
NNN is a series number that increases with every executable xml block in the input file.
\subsection{Time-step Study}
Now we will study the dependence of the DMC energy with time-step. From the top directory,
go to “ex1\_first-run-hartree-fock/dmc\_timestep”. This folder contains a basic xml input
file (dmc\_ts.xml) that performs a short VMC calculation and three DMC calculations
with varying time-steps (0.1, 0.05, 0.01). Link the particle set and the last optimization
file from the previous folder (the file called jopt-h2o.sNNN.opt.xml with the largest value of
NNN). Rename the optimized wave-function to any suitable name if you wish, for example
h2o.opt.xml, and change the name of the particle set and wave-function files in the
input file. An optimized wave-function can be found in the reference files (same location)
\subsection{Time-step study}
Now we will study the dependence of the DMC energy with time step. From the top directory,
go to “\texttt{ex1\_first-run-hartree-fock/dmc\_timestep}.” This folder contains a basic XML input
file (\texttt{dmc\_ts.xml}) that performs a short VMC calculation and three DMC calculations
with varying time steps (0.1, 0.05, 0.01). Link the \texttt{particleset} and the last \texttt{optimization}
file from the previous folder (the file called \texttt{jopt-h2o.sNNN.opt.xml} with the largest value of
NNN). Rename the optimized \texttt{wavefunction} file to any suitable name if you wish (for example,
\texttt{h2o.opt.xml}) and change the name of the \texttt{particleset} and \texttt{wavefunction} files in the
input file. An optimized wavefunction can be found in the reference files (same location)
in case it is needed. %Using the submission script of the previous exercise as a base, create a
%submission script for this step and submit the run. Set the number of nodes to 32 (2 places
%must be changed), the number of threads to 16 and leave the number of tasks at 1.
@ -202,14 +202,14 @@ cp ../opt/jopt-h2o.s007.opt.xml h2o.opt.wfs.xml
# edit dmc_ts.xml to include the correct ptcl.xml and wfs.xml
jobrun_vesta qmcpack dmc_ts.xml
\end{lstlisting}
While these runs complete, go to section D and review the basic VMC and DMC input
blocks. Notice that in the current DMC blocks, as the time-step is decreased the number of blocks is also increased. Why is this?
While these runs complete, go to Section\ref{sec:lab_adv_mol_vmcdmc_appendix}, ``Appendix D: VMC and DMC XML block," and review the basic VMC and DMC input
blocks. Notice that in the current DMC blocks the time step is decreased as the number of blocks is increased. Why is this?
When the simulations are finished, use qmca to analyze the output files and to plot the
DMC energy as a function of time-step. Results should be qualitatively similar to those
presented in figure~\ref{fig:lam_dmc_timestep}, in this case we present more time-steps with well converged results to
better illustrate the time-step dependence. In realistic calculations, the time-step must be
chosen small enough so that the resulting error is below the desire accuracy. Alternatively,
When the simulations are finished, use \texttt{qmca} to analyze the output files and plot the
DMC energy as a function of time step. Results should be qualitatively similar to those
presented in Figure \ref{fig:lam_dmc_timestep}; in this case we present more time steps with well converged results to
better illustrate the time step dependence. In realistic calculations, the time step must be
chosen small enough so that the resulting error is below the desired accuracy. Alternatively,
various calculations can be performed and the results extrapolated to the zero time-step
limit.
@ -218,35 +218,35 @@ limit.
\begin{center}
\includegraphics[trim = 0mm 0mm 0mm 0mm, clip,width=0.75\columnwidth]{./figures/lab_advanced_molecules_dmc_timestep.png}
\end{center}
\caption{DMC energy as a function of timestep.}
\caption{DMC energy as a function of time step.}
\label{fig:lam_dmc_timestep}
\end{figure}
\subsection{Walker Population Study}
\subsection{Walker population study}
Now we will study the dependence of the DMC energy with the number of walkers in the
simulation. Remember that, in principle, the DMC distribution is reached in the limit of
an infinite number of walkers. In practice, the energy and most properties converge to high
accuracy with $\sim$100-1000 walkers. The actual number of walkers needed in a calculation
will depend on the accuracy of the VMC wave-function and on the complexity and size of
the system. Also notice that using too many walkers is not a problem, at worse it will be
accuracy with $\sim$100--1,000 walkers. The actual number of walkers needed in a calculation
will depend on the accuracy of the VMC wavefunction and on the complexity and size of
the system. Also notice that using too many walkers is not a problem; at worse it will be
inefficient since it will cost more computer time than necessary. In fact, this is the strategy
used when running QMC calculations on large parallel computers since we can reduce the
statistical error bars efficiently by running with large walker populations distributed across
all processors.
From the top directory, go to ``ex1\_first-run-hartree-fock/dmc\_walkers''. Copy the
optimized wave-function and particle set files used in the previous calculations to the current
folder, these are the ones generated on step 2 of this exercise. An optimized wave-function can be found in the reference files (same location) in case it is needed. The directory
contains a sample DMC input file and submission script. Make 3 directories named NWx,
with x values 120,240,480 and copy the input file to each one. Go
to ``NW120'', and, in the input file, change the name of the wave-function and particle set
files (in this case they will be located one directory above, so use ``../dmc\_timestep/h2.opt.xml'' for
example), change the pseudopotential directory to point to one directory above, change ``targetWalkers'' to 120, change the number of steps to 100, the time-step
to 0.01 and the number of blocks to 400. Notice that ``targetWalkers'' is one way to set the desired (average) number of walkers in a DMC calculation. One can alternatively set ``samples'' in the \ixml{<qmc method=``vmc''} block to carry over de-correlated VMC configurations as DMC walkers. For your own simulations we generally recommend setting $\sim$2*(\#threads)
From the top directory, go to ``\texttt{ex1\_first-run-hartree-fock/dmc\_walkers}.'' Copy the
optimized \texttt{wavefunction} and \texttt{particleset} files used in the previous calculations to the current
folder; these are the files generated during step 2 of this exercise. An optimized \texttt{wavefunction} file can be found in the reference files (same location) in case it is needed. The directory
contains a sample DMC input file and submission script. Create three directories named NWx,
with x values of 120,240,480, and copy the input file to each one. Go
to ``NW120,'' and, in the input file, change the name of the \texttt{wavefunction} and \texttt{particleset}
files (in this case they will be located one directory above, so use ``\texttt{../dmc\_timestep/h2.opt.xml},'' for
example); change the PP directory so that it points to one directory above; change ``targetWalkers'' to 120; and change the number of steps to 100, the time step
to 0.01, and the number of blocks to 400. Notice that ``targetWalkers'' is one way to set the desired (average) number of walkers in a DMC calculation. One can alternatively set ``samples'' in the \ixml{<qmc method="vmc"} block to carry over de-correlated VMC configurations as DMC walkers. For your own simulations, we generally recommend setting $\sim$2*(\#threads)
walkers per node (slightly smaller than this value).
The main steps needed to perform this exercise are:
The main steps needed to perform this exercise are
\begin{shade}
cd ${TRAINING TOP}/ex1_first-run-hartree-fock/dmc_walkers
cp ../opt/h2o.ptcl.xml ./
@ -272,58 +272,58 @@ jobrun_vesta qmcpack dmc_wk.xml
Repeat the same procedure in the other folders by setting (targetWalkers=240,
steps=100, timestep=0.01, blocks=200) in NW240 and (targetWalkers=480,
steps=100, timestep=0.01, blocks=100) in NW480. When
the simulations complete, use qmca to analyze and plot the energy as a function of the
number of walkers in the calculation. As always, figure~\ref{fig:lam_dmc_popcont}
the simulations complete, use \texttt{qmca} to analyze and plot the energy as a function of the
number of walkers in the calculation. As always, Figure \ref{fig:lam_dmc_popcont}
shows representative results of the
dependence of the energy on the number of walkers for a single water molecule. As shown,
energy dependence on the number of walkers for a single water molecule. As shown,
less than 240 walkers are needed to obtain an accuracy of 0.1 mHa.
\section{Exercise \#2 Slater-Jastrow Wave-Function Options}
\section{Exercise \#2 Slater-Jastrow wavefunction options}
From this point on in the tutorial we assume familiarity with the basic parameters in the
optimization, VMC and DMC XML input blocks of QMCPACK. In addition, we assume
optimization, VMC, and DMC XML input blocks of QMCPACK. In addition, we assume
familiarity with the submission system. As a result, the folder structure will not contain
any prepared input or submission files, the student will generate them using
any prepared input or submission files, so you will need to generate them using
input files from exercise 1. In the case of QMCPACK sample
files, you will find optm.xml, vmc dmc.xml and submit.csh files. Some of
files, you will find \texttt{optm.xml}, \texttt{vmc dmc.xml}, and \texttt{submit.csh files}. Some of
the options in these files can be left unaltered, but many of them will need to be tailored to
the particular calculation.
In this exercise we will study the dependence of the DMC energy on the choices made
in the wave-function ansatz. In particular, we will study the influence/dependence of the
in the wavefunction ansatz. In particular, we will study the influence/dependence of the
VMC energy with the various terms in the Jastrow. We will also study the influence of
the VMC and DMC energies on the single particle orbitals used to form the Slater determinant
in single determinant wave-functions. For this we will use wave-functions generated
the VMC and DMC energies on the SPOs used to form the Slater determinant
in single-determinant wavefunctions. For this we will use wavefunctions generated
with various exchange-correlation functionals in DFT. Finally, we will optimize a simple
multi-determinant wave-function and study the dependence of the energy o the number of
multideterminant wavefunction and study the dependence of the energy on the number of
configurations used in the expansion. All of these exercises will be performed on the water
molecule at equilibrium.
\subsection{Influence of Jastrow on VMC energy with HF wave-function}
\subsection{Influence of Jastrow on VMC energy with HF wavefunction}
In this section we will study the dependence of the VMC energy on the various Jastrow
terms, e.g. one-body, two-body and three-body. From the top directory, go to ``ex2\_slater-jastrow-wf-options/jastrow''.
We will compare the single determinant VMC energy using a two-body
Jastrow term, both one- and two-body terms and finally one-, two- and three-body
terms (e.g., 1-body, 2-body and 3-body. From the top directory, go to ``\texttt{ex2\_slater-jastrow-wf-options/jastrow.''}
We will compare the single-determinant VMC energy using a 2-body
Jastrow term, both 1- and 2-body terms, and finally 1-, 2- and 3-body
terms. Since we are interested in the influence of the Jastrow, we will use the HF orbitals
calculated in exercise \#1. Make three folders named 2j,12j,123j. For both 2j and
12j %(we have already optimized a wave-function for the 1-2-3J case, so the steps will be
calculated in exercise \#1. Make three folders named 2j, 12j, and 123j. For both 2j and
12j, %(we have already optimized a wave-function for the 1-2-3J case, so the steps will be
%slightly different in this case)
, copy the input file optm.xml %and the sample submission file
from ``ex1\_first-run-hartree-fock/opt'' . This input file performs both wave-function optimization
and a VMC calculation. Remember to correct relative paths to the pseudopotential directory. Copy the un-optimized HF wave-function and particle set files
from ``ex1\_first-run-hartree-fock/convert'', if you followed the instructions in exercise \#1 these should be
named h2o.wfs.xml and h2o.ptcl.xml. Otherwise, you can obtained them from the
REFERENCE files. Modify the file h2o.wfs.xml to remove the appropriate jastrow
blocks. For example, for a two-body Jastrow (only), you need to eliminate the jastrow
blocks named \ixml{<jastrow name="J1"} and \ixml{<jastrow name="J3"}. In the case of 12j, remove
only \ixml{<jastrow name="J3"}. Recommended settings for the optimization run are: nodes=32,
threads=16, blocks=250, samples=128000, time-step=0.5, 8 optimization loops, and in the
VMC section we recommend walkers=16, blocks=1000, steps=1, substeps=100. Notice that
copy the input file \texttt{optm.xml} %and the sample submission file
from ``\texttt{ex1\_first-run-hartree-fock/opt.}'' This input file performs both wavefunction optimization
and a VMC calculation. Remember to correct relative paths to the PP directory. Copy the un-optimized HF \texttt{wavefunction} and \texttt{particleset} files
from ``\texttt{ex1\_first-run-hartree-fock/convert}''; if you followed the instructions in exercise \#1 these should be
named \texttt{h2o.wfs.xml} and \texttt{h2o.ptcl.xml}. Otherwise, you can obtained them from the
REFERENCE files. Modify the \texttt{h2o.wfs.xml} file to remove the appropriate Jastrow
blocks. For example, for a 2-body Jastrow (only), you need to eliminate the Jastrow
blocks named \ixml{<jastrow name="J1"} and \ixml{<jastrow name="J3."} In the case of 12j, remove
only \ixml{<jastrow name="J3."} Recommended settings for the optimization run are nodes=32,
threads=16, blocks=250, samples=128000, time-step=0.5, 8 optimization loops. Recommended settings in the
VMC section are walkers=16, blocks=1000, steps=1, substeps=100. Notice that
samples should always be set to blocks*threads per node*nodes = 32*16*250=128000. Repeat
the process in both 2j and 12j cases. For the 123j case, the wave-function has
already been optimized in the previous exercise. Copy the optimized HF wave-function and
the particle set from ``ex1\_first-run-hartree-fock/opt''. Copy the input file from any of the previous runs and remove the optimization block from the
the process in both 2j and 12j cases. For the 123j case, the wavefunction has
already been optimized in the previous exercise. Copy the optimized HF wavefunction and
the particleset from ``\texttt{ex1\_first-run-hartree-fock/opt.}'' Copy the input file from any of the previous runs and remove the optimization block from the
input, just leave the VMC step. In all three cases, modify the submission script and submit the run.
\begin{figure}
@ -334,71 +334,70 @@ input, just leave the VMC step. In all three cases, modify the submission script
\label{fig:lam_vmc_jastrow}
\end{figure}
These simulations will take several minutes to complete. This is an excellent opportunity
to go to section E and review the wavefunction XML block used by QMCPACK. When the
simulation are completed, use qmca to analyze the output files. Using your favorite plotting
program (e.g. gnu plot), plot the energy and variance as a function of the Jastrow form.
Figure~\ref{fig:lam_vmc_jastrow} shows a typical result for this calculation. As can be seen, the VMC energy and
Because these simulations will take several minutes to complete, this is an excellent opportunity
to go to Section~\ref{sec:lab_adv_mol_wf_appendix}, ``Appendix E: Wavefunction XML block," and review the wavefunction XML block used by QMCPACK. When the
simulations are completed, use \texttt{qmca} to analyze the output files. Using your favorite plotting
program (e.g., gnu plot), plot the energy and variance as a function of the Jastrow form.
Figure \ref{fig:lam_vmc_jastrow} shows a typical result for this calculation. As can be seen, the VMC energy and
variance depends strongly on the form of the Jastrow. Since the DMC error bar is directly
related to the variance of the VMC energy, improving the Jastrow will always lead to a
reduction in the DMC effort. In addition, systematic approximations (time-step, number of
walkers, etc) are also reduced with improved wave-functions.
reduction in the DMC effort. In addition, systematic approximations (time step, number of
walkers, etc.) are also reduced with improved wavefunctions.
\subsection{Generation of wave-functions from DFT using GAMESS}
In this section we will use GAMESS to generate wave-functions for QMCPACK from
DFT calculations. From the top folder, go to ``ex2\_slater-jastrow-wf-options/orbitals''. In order to demonstrate
\subsection{Generation of wavefunctions from DFT using GAMESS}
In this section we will use GAMESS to generate wavefunctions for QMCPACK from
DFT calculations. From the top folder, go to ``\texttt{ex2\_slater-jastrow-wf-options/orbitals}.'' To demonstrate
the variation in DMC energies with the choice of DFT orbitals, we will choose the following
set of exchange-correlation functionals (PBE, PBE0, BLYP, B3LYP). For each functional,
make a directory using your preferred naming convention (e.g. the name of the functional).
make a directory using your preferred naming convention (e.g., the name of the functional).
Go into each folder and copy a GAMESS input file from %for a ROHF calculation from
``ex1\_first-run-hartree-fock/gms'' .%, a file named rohf.inp should exist.
Rename the file with your preferred naming convention, we suggest using h2o.[dft].inp, where [dft] is the name of
``\texttt{ex1\_first-run-hartree-fock/gms}.'' %, a file named rohf.inp should exist.
Rename the file with your preferred naming convention; we suggest using \texttt{h2o.[dft].inp}, where [dft] is the name of
the functional used in the calculation. At this point, this input file should be identical to the
one used to generate the HF wave-function in exercise \#1. In order to perform a DFT
one used to generate the HF wavefunction in exercise \#1. To perform a DFT
calculation we only need to add ``DFTTYP'' to the \igamess{\$CONTRL ... \$END} section and set
it to the desired functional type, for example ``DFTTYP=PBE'' for a PBE functional. This
it to the desired functional type, for example, ``DFTTYP=PBE'' for a PBE functional. This
variable must be set to (PBE, PBE0, BLYP, B3LYP) to obtain the appropriate functional in
GAMESS. For a complete list of implemented functionals, see the GAMESS input manual.
\subsection{Optimization and DMC calculations with DFT wave-functions}
In this section we will optimize the wave-function generated in the previous step and
perform DMC calculations. From the top directory, go to “ex2\_slater-jastrow-wf-options/orbitals”.
The steps required to achieve this are identical to those used to optimize the wave-function
\subsection{Optimization and DMC calculations with DFT wavefunctions}
In this section we will optimize the wavefunction generated in the previous step and
perform DMC calculations. From the top directory, go to “\texttt{ex2\_slater-jastrow-wf-options/orbitals}.”
The steps required to achieve this are identical to those used to optimize the wavefunction
with HF orbitals. Make individual folders for each calculation and obtain the necessary files
to perform optimization, VMC and DMC calculations from ``ex1\_first-run-hartree-fock/opt'' and ``ex1\_first-run-hartree-fock/dmc\_ts'', for example.
to perform optimization, for example, VMC and DMC calculations from ``for \texttt{ex1\_first-run-hartree-fock/opt}'' and ``\texttt{ex1\_first-run-hartree-fock/dmc\_ts}.''
%A file named optm vmc dmc.xml should exist that contains all three execution blocks.
For each functional, make the appropriate modifications to the input files and copy the particle
set and wave-function files from the appropriate directory in “ex2\_slater-jastrow-wf-options/orbitals/[dft]”. We
For each functional, make the appropriate modifications to the input files and copy the \texttt{particleset} and \texttt{wavefunction} files from the appropriate directory in “\texttt{ex2\_slater-jastrow-wf-options/orbitals/[dft]}.” We
recommend the following settings: nodes=32, threads=16, (in optimization) blocks=250,
samples=128000, timestep=0.5, 8 optimization loops, (in VMC) walkers=16, blocks=100,
steps=1, substeps=100, (in DMC) blocks 400, targetWalkers=960, timestep=0.01. Submit
the runs and analyze the results using qmca .
steps=1, substeps=100, (in DMC) blocks 400, targetWalkers=960, and timestep=0.01. Submit
the runs and analyze the results using \texttt{qmca}.
How do the energies compare against each other? How do they compare against DMC
energies with HF orbitals?
%Orbital Sets and Configurations in
\section{Exercise \#3: Multi-Determinant Wave-Functions}
\section{Exercise \#3: Multideterminant wavefunctions}
In this exercise we will study the dependence of the DMC energy on the set of orbitals
and the type of configurations included in a multi-determinant wave-function.
and the type of configurations included in a multideterminant wavefunction.
\subsection{Generation of a CISD wave-functions using GAMESS}
In this section we will use GAMESS to generate a multi-determinant wave-function with
Configuration Interaction with Single and Double excitations (CISD). In CISD, the Schrodinger equation is solved exactly in a basis of determinants
\subsection{Generation of a CISD wavefunctions using GAMESS}
In this section we will use GAMESS to generate a multideterminant wavefunction with
configuration interaction with single and double excitations (CISD). In CISD, the Schrodinger equation is solved exactly on a basis of determinants
including the HF determinant and all its single and double excitations.
Go to ``ex3\_multi-slater-jastrow/cisd/gms'' and you'll see input and output files named h2o.cisd.inp and h2o.cisd.out. Due to technical problems with GAMESS in the BGQ architecture of VESTA, we are unable to use CISD properly in GAMESS. For this reason, the output of the calculation is already provided in the directory.
Go to ``\texttt{ex3\_multi-slater-jastrow/cisd/gms}'' and you will see input and output files named \texttt{h2o.cisd.inp} and \texttt{h2o.cisd.out}. Because of technical problems with GAMESS in the BGQ architecture of VESTA, we are unable to use CISD properly in GAMESS. Consequently, the output of the calculation is already provided in the directory.
%You'll see several input and output files named h2o.XXX.inp
%and h2o.XXX.out, where XXX is one of the following multi-determinant methods: CISD,
%CASSCF, CASCI, SOCI.
There will be time in the next step to study the GAMESS input
files and the description in section A. %In the next exercise we will use the CISD output, in
files and the description in Section~\ref{sec:lab_adv_mol_gamess}, ``Appendix A: GAMESS input." %In the next exercise we will use the CISD output, in
%the next exercise we will use the remaining files.
Since the output is already provided, the
only thing needed is to use the converter to generate the appropriate QMCPACK files. %Copy a submission script from GAMESS/Generic Files and execute the converter for all the output
only action needed is to use the converter to generate the appropriate QMCPACK files. %Copy a submission script from GAMESS/Generic Files and execute the converter for all the output
%files in the directory (with the exception of CASSCF, which is used to generate orbitals).
%but it doesn’t contain appropriate CI coefficients).
\begin{shade}
@ -407,27 +406,27 @@ jobrun_vesta convert4qmc h2o.cisd.out -ci h2o.cisd.out \
\end{shade}
We used the PRTMO=.T. flag in the GUESS section to include orbitals in the output file. You should read these orbitals from the output (-readInitialGuess 40).
The highest occupied orbital in any determinant should be 34, so reading 40 orbitals is a safe choice. In this case, it is important to rename the xml files with meaningful names, for example h2o.cisd.wfs.xml. A threshold of 0.0075 is sufficient for the calculations in the training.
The highest occupied orbital in any determinant should be 34, so reading 40 orbitals is a safe choice. In this case, it is important to rename the XML files with meaningful names, for example, \texttt{h2o.cisd.wfs.xml}. A threshold of 0.0075 is sufficient for the calculations in the training.
\subsection{Optimization of Multi-Determinant wave-function}
\subsection{Optimization of a multideterminant wavefunction}
In this section we will optimize the wave-function generated in the previous step. There
is no difference in the optimization steps if a single determinant and a multi-determinant wave-function.
QMCPACK will recognize the presence of a multi-determinant wavefunction and will automatically
optimize the linear coefficients by default. Go to ``ex3\_multi-slater-jastrow/cisd'' and make a folder called
thres0.01. Copy the particle set and wavefunction files created in the previous step to the current
directory. With your favorite text editor, open the wave-function file h2o.wfs.xml. Look for
In this section we will optimize the wavefunction generated in the previous step. There
is no difference in the optimization steps if a single determinant and a multideterminant wavefunction.
QMCPACK will recognize the presence of a multideterminant wavefunction and will automatically
optimize the linear coefficients by default. Go to ``\texttt{ex3\_multi-slater-jastrow/cisd}'' and make a folder called
\texttt{thres0.01}. Copy the \texttt{particleset} and \texttt{wavefunction} files created in the previous step to the current
directory. With your favorite text editor, open the \texttt{wavefunction} file \texttt{h2o.wfs.xml}. Look for
the multideterminant XML block and change the ``cutoff'' parameter in detlist to 0.01. Then follow
the same steps used in the subsection ``Optimization and DMC calculations with DFT wave-functions''
to optimize the wave-function. Similar to this case, design a QMCPACK input file that performs
wave-function optimization followed by VMC and DMC calculations. Submit the calculation.
the same steps used in Section 9.4.3, ``Optimization and DMC calculations with DFT wavefunctions''
to optimize the wavefunction. Similar to this case, design a QMCPACK input file that performs
wavefunction optimization followed by VMC and DMC calculations. Submit the calculation.
This is a good time to review the GAMESS input file description in Appendix section A.
This is a good time to review the GAMESS input file description in Section~\ref{sec:lab_adv_mol_gamess}, ``Appendix A. GAMESS input."
When the run is completed, go to the previous directory and make a new folder named
thres0.0075. Repeat the steps performed above to optimize the wave-function with a cutoff of 0.01, but use a cutoff of 0.0075 this time. This will increase the number of determinants used in the calculation. Notice the ``cutoff'' parameter in the XML should be less than the ``-threshold 0.0075'' flag passed to the converted, which is further bounded by the PRTTOL flag in the GAMESS input.
\texttt{thres0.0075}. Repeat the previous steps to optimize the wavefunction with a cutoff of 0.01, but use a cutoff of 0.0075 this time. This will increase the number of determinants used in the calculation. Notice that the ``cutoff'' parameter in the XML should be less than the ``-threshold 0.0075'' flag passed to the converted, which is further bounded by the PRTTOL flag in the GAMESS input.
After the wave-function is generated, we are ready to optimize. Instead of starting from an un-optimized wave-function, we can start from optimized wave-function from thres0.01 to speed up convergence. You will need to modify the file and change the cutoff in detlist to 0.0075 with a text editor. Repeat the optimization steps and submit the calculation.
After the wavefunction is generated, we are ready to optimize. Instead of starting from an un-optimized wavefunction, we can start from the optimized wavefunction from thres0.01 to speed up convergence. You will need to modify the file and change the cutoff in detlist to 0.0075 with a text editor. Repeat the optimization steps and submit the calculation.
\begin{figure}
\begin{center}
@ -437,56 +436,54 @@ After the wave-function is generated, we are ready to optimize. Instead of start
\label{fig:lam_dmc_ci_cisd}
\end{figure}
When you are done, use qmca to analyze the results. Compare the energies at these two
coefficient cutoffs with the energies obtained with DFT orbitals. Due to the time limitations of this tutorial it is not practical to optimize the wave-functions with a smaller cutoff, since this would require more samples and longer runs due to the larger number of optimizable parameters. Figure~\ref{fig:lam_dmc_ci_cisd} shows the results of such exercise, the DMC energy as a function of the cutoff in the wave-function. As can be seen, a large improvement in the energy is obtained as the number of configurations is increased.
When you are done, use \texttt{qmca} to analyze the results. Compare the energies at these two
coefficient cutoffs with the energies obtained with DFT orbitals. Because of the time limitations of this tutorial, it is not practical to optimize the wavefunctions with a smaller cutoff since this would require more samples and longer runs due to the larger number of optimizable parameters. Figure \ref{fig:lam_dmc_ci_cisd} shows the results of such exercise: the DMC energy as a function of the cutoff in the wavefunction. As can be seen, a large improvement in the energy is obtained as the number of configurations is increased.
%Since the un-optimized wave-functions were generated in subsection “Generation of a CISD wavefunctions
%using GAMESS” of exercise \#2, we can skip this section and go straight to the
%wave-function optimization.
\subsection{CISD, CASCI and SOCI}
\subsection{CISD, CASCI, and SOCI}
Go to “ex3\_multi-slater-jastrow” and inspect folders for the remaining wave-function types: CASCI and SOCI. Follow steps in the previous exercise and obtain the optimized wave-functions for these determinant choices. Notice the SOCI GAMESS output is not included because it is large. Already converted XML inputs can be found in ``ex3\_multi-slater-jastrow/soci/thres*''. %The exercise has already been performed with a CISD wave-function in exercise \#2.
Go to “\texttt{ex3\_multi-slater-jastrow}” and inspect the folders for the remaining wavefunction types: CASCI and SOCI. Follow the steps in the previous exercise and obtain the optimized wavefunctions for these determinant choices. Notice that the SOCI GAMESS output is not included because it is large. Already converted XML inputs can be found in ``\texttt{ex3\_multi-slater-jastrow/soci/thres*}.'' %The exercise has already been performed with a CISD wave-function in exercise \#2.
A CASCI wave-function is produced from a CI calculation that includes all the determinants
A CASCI wavefunction is produced from a CI calculation that includes all the determinants
in a complete active space (CAS) calculation, in this case using the orbitals from a previous CASSCF
calculation. In this case we used a CAS(8,8) active space, that includes all determinants
generated by distributing 8 electrons in the lowest 8 orbitals. A second-order CI (SOCI) calculation is similar
calculation. In this case we used a CAS(8,8) active space that includes all determinants
generated by distributing 8 electrons in the lowest 8 orbitals. A SOCI calculation is similar
to the CAS-CI calculation, but in addition to the determinants in the CAS it also includes
all single and double excitations from all of them, leading to a much larger determinant
set. Since we now have considerable experience optimizing wave-functions and calculating
DMC energies, we will leave it to the student to complete the remaining tasks on its own.
set. Since you now have considerable experience optimizing wavefunctions and calculating
DMC energies, we will leave it to you to complete the remaining tasks on your own.
If you need help, refer to previous exercises in the tutorial. Perform optimizations for both
wave-functions using cutoffs in the CI expansion of 0.01 an 0.0075. If there is enough time
left, try to optimize the wave-functions with a cutoff of 0.005. Analyze the results and plot
the energy as a function of cutoff for all three cases, CISD, CAS-CI and SOCI.
wavefunctions using cutoffs in the CI expansion of 0.01 an 0.0075. If you have time, try to optimize the wavefunctions with a cutoff of 0.005. Analyze the results and plot
the energy as a function of cutoff for all three cases: CISD, CAS-CI, and SOCI.
Figure \ref{fig:lam_dmc_ci_cisd} shows the result of similar calculations using more samples and smaller cutoffs.
The results should be similar to those produced in the tutorial. For reference, the exact
energy of the water molecule with ECPs is approximately -17.276 Ha. From the results of the
tutorial, how does the selection of determinants is related to the expected DMC energy?
tutorial, how does the selection of determinants relate to the expected DMC energy?
What about the choice in the set of orbitals?
\newpage
\section{Appendix A: GAMESS input}
\section{Appendix A: GAMESS input}\label{sec:lab_adv_mol_gamess}
In this section we provide a brief description of the GAMESS input needed to produce
trial wave-function for QMC calculations with QMCPACK. We assume basic familiarity
with GAMESS input structure, in particular regarding the input of atomic coordinates and
the definition of gaussian basis sets. This section will focus on the generation of the output
files needed by the converter tool, convert4qmc. For a description of the converter, see B.
trial wavefunction for QMC calculations with QMCPACK. We assume basic familiarity
with GAMESS input structure, particularly regarding the input of atomic coordinates and
the definition of Gaussian basis sets. This section focuses on generation of the output
files needed by the converter tool, \texttt{convert4qmc}. For a description of the converter, see Section~\ref{sec:lab_adv_mol_convert4qmc}, ``Appendix B: convert4qmc."
Only a subset of the methods available in GAMESS can be used to generate wave-functions
for QMCPACK and we restrict our description here to these.
Only a subset of the methods available in GAMESS can be used to generate wavefunctions
for QMCPACK, and we restrict our description to these.
For a complete description of all the options and methods available
in GAMESS, please refer to the official documentation which could be found in
”http://www.msg.ameslab.gov/gamess/documentation.html”.
in GAMESS, please refer to the official documentation at ``http://www.msg.ameslab.gov/gamess/documentation.html.”
Currently, convert4qmc can process output for the following methods in GAMESS (in
SCFTYP) : RHF, ROHF, and MCSCF. Both HF as well as DFT calculations (any DFT
type) could be used in combination with RHF and ROHF calculations. For MCSCF and CI
calculations, ALDET, ORMAS and GUGA drivers can be used (see below for details).
Currently, \texttt{convert4qmc} can process output for the following methods in GAMESS (in
SCFTYP): RHF, ROHF, and MCSCF. Both HF and DFT calculations (any DFT
type) can be used in combination with RHF and ROHF calculations. For MCSCF and CI
calculations, ALDET, ORMAS, and GUGA drivers can be used (details follow).
\subsection{HF input}
@ -510,9 +507,9 @@ $END
Main options:
\begin{enumerate}
\item{SCFTYP: Type of SCF method, options: RHF, ROHF, MCSCF, UHF and NONE.}
\item{RUNTYP: Type of run. For QMCPACK wave-function generation this should always be ENERGY.}
\item{RUNTYP: Type of run. For QMCPACK wavefunction generation this should always be ENERGY.}
\item{MULT: Multiplicity of the molecule.}
\item{ISPHER: Use spherical harmonics (1) or cartesian basis functions (-1).}
\item{ISPHER: Use spherical harmonics (1) or Cartesian basis functions (-1).}
\item{COORD: Input structure for the atomic coordinates in \$DATA.}
\end{enumerate}
@ -521,17 +518,17 @@ Main options:
The main difference between the input for a RHF/ROHF calculation and a DFT calculation
is the definition of the DFTTYP parameter. If this is set in the \$CONTROL
section, a DFT calculation will be performed with the appropriate functional. Notice that
while the default values are usually adequate, DFT calculations have many options involving
although the default values are usually adequate, DFT calculations have many options involving
the integration grids and accuracy settings. Make sure you study the input manual to be
aware of these. Refer to the input manual for a list of the implemented exchange-correlation
functionals.
\subsection{Multi-Configuration Self-Consistent Field (MCSCF)}
MCSCF calculations are performed by setting SCFTYP=MCSCF in the $CONTROL
section. If this option is set, a $MCSCF section must be added to the input file with the
\subsection{MCSCF}
MCSCF calculations are performed by setting SCFTYP=MCSCF in the CONTROL
section. If this option is set, an MCSCF section must be added to the input file with the
options for the calculation. An example section for the water molecule used in the tutorial
is shown below.
follows.
\begin{lstlisting}[style=GAMESS]
$MCSCF CISTEP=GUGA MAXIT=1000 FULLNR=.TRUE. ACURCY=1.0D-5 $END
@ -542,22 +539,21 @@ options compatible with QMCPACK are: ALDET, GUGA, and ORMAS. Depending on the
package used, additional input sections are needed.
\subsection{Configuration Interaction (CI)}
\subsection{CI}
Configuration interaction (full CI, truncated CI, CAS-CI, etc) calculations are performed
by setting \igamess{SCFTYP=NONE} and \igamess{CITYP=GUGA,ALDET,ORMAS}. Each one of this packages
requires further input sections, which are typically slightly different to the input sections
by setting \igamess{SCFTYP=NONE} and \igamess{CITYP=GUGA,ALDET,ORMAS}. Each one of these packages
requires further input sections, which are typically slightly different from the input sections
needed for MCSCF runs.
\subsection{GUGA: Unitary Group CI package}
The GUGA package is the only alternative if one wants CSFs with GAMESS. Below
we provide a very brief description of input sections needed to perform MCSCF, CASCI,
truncated CI and SOCI with this package. For a complete description of these methods and
\subsection{GUGA: Unitary group CI package}
The GUGA package is the only alternative if one wants CSFs with GAMESS. We subsequently provide a very brief description of the input sections needed to perform MCSCF, CASCI,
truncated CI, and SOCI with this package. For a complete description of these methods and
all the options available, please refer to the GAMESS input manual.
\subsubsection{GUGA-MCSCF}
The following input section performs a CASCI calculation, with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), e.g. CAS(8,8). NMCC is the number of frozen
The following input section performs a CASCI calculation with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), for example, CAS(8,8). NMCC is the number of frozen
orbitals (doubly occupied orbitals in all determinants), NDOC is the number of double
occupied orbitals in the reference determinant, NVAL is the number of singly occupied
orbitals in the reference (for spin polarized cases), and NVAL is the number of orbitals in
@ -570,8 +566,8 @@ $DRT GROUP=C2v NMCC=0 NDOC=4 NALP=0 NVAL=4 ISTSYM=1 MXNINT= 500000 FORS=.TRUE. $
\end{lstlisting}
\subsubsection{GUGA-CASCI}
The following input section performs a CASCI calculation, with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), e.g. CAS(8,8). NFZC is the number of frozen
The following input section performs a CASCI calculation with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), for example, CAS(8,8). NFZC is the number of frozen
orbitals (doubly occupied orbitals in all determinants). All other parameters are identical
to those in the MCSCF input section.
@ -580,10 +576,10 @@ $CIDRT GROUP=C2v NFZC=0 NDOC=4 NALP=0 NVAL=4 NPRT=2 ISTSYM=1 FORS=.TRUE. MXNINT=
$GUGDIA PRTTOL=0.001 CVGTOL=1.0E-5 ITERMX=1000 $END
\end{lstlisting}
\subsubsection{GUGA-Truncated CI}
The following input sections will lead to a truncated CI calculation, in this particular case
\subsubsection{GUGA-truncated CI}
The following input sections will lead to a truncated CI calculation. In this particular case
it will perform a CISD calculation since IEXCIT is set to 2. Other values in IEXCIT will lead
to different CI truncations, for example IEXCIT=4 will lead to CISDTQ. Notice that only
to different CI truncations; for example, IEXCIT=4 will lead to CISDTQ. Notice that only
the lowest 30 orbitals will be included in the generation of the excited determinants in this
case. For a full CISD calculation, NVAL should be set to the total number of virtual orbitals.
@ -593,8 +589,8 @@ $GUGDIA PRTTOL=0.001 CVGTOL=1.0E-5 ITERMX=1000 $END
\end{lstlisting}
\subsubsection{GUGA-SOCI}
The following input section performs a SOCI calculation, with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), e.g. CAS(8,8). Since SOCI is set to .TRUE.,
The following input section performs a SOCI calculation with a CAS that includes 8
electrons in 8 orbitals (4 DOC and 4 VAL), for example, CAS(8,8). Since SOCI is set to .TRUE.,
all single and double determinants from all determinants in the CAS(8,8) will be included.
\begin{lstlisting}[style=GAMESS]
@ -604,10 +600,10 @@ $GUGDIA PRTTOL=0.001 CVGTOL=1.0E-5 ITERMX=1000 $END
\subsection{ECP}
To use Effective Core Potentials (ECP) in GAMESS, you must define a \{\ishell{\$ECP ... \$END}\}
To use ECPs in GAMESS, you must define a \{\ishell{\$ECP ... \$END}\}
block. There must be a definition of a potential for every atom in the system, including
symmetry equivalent ones. In addition, they must appear in the particular order expected
by GAMESS. Below is an example of an ECP input block for a single water molecule using
by GAMESS. The following example shows an ECP input block for a single water molecule using
BFD ECPs. To turn on the use of ECPs, the option “ECP=READ” must be added to the
CONTROL input block.
@ -631,26 +627,26 @@ $END
\newpage
\section{Appendix B: convert4qmc}
\section{Appendix B: convert4qmc}\label{sec:lab_adv_mol_convert4qmc}
To generate the particleset and wavefunction XML blocks required by QMCPACK in
calculations with molecular systems, the converter convert4qmc must be used. The converter
will read the standard output from the appropriate Quantum Chemistry calculation and will
generate all the necessary input for QMCPACK. Below we describe the main options of the
converter for GAMESS output. In general, there are 3 ways to use the converter depending
on the type of calculation performed. The minimum syntax for each option is found below.
For a description of the xml files produced by the converter, see section E.
calculations with molecular systems, the converter \texttt{convert4qmc} must be used. The converter
will read the standard output from the appropriate quantum chemistry calculation and will
generate all the necessary input for QMCPACK. In the following, we describe the main options of the
converter for GAMESS output. In general, there are three ways to use the converter depending
on the type of calculation performed. The minimum syntax for each option is shown subsequently.
For a description of the XML files produced by the converter, see Section~\ref{sec:lab_adv_mol_wf_appendix}, ``Appendix E: Wavefunction XML block."
\begin{enumerate}
\item{For all single determinant calculations (HF and DFT with any DFTTYP):}
\item{For all single-determinant calculations (HF and DFT with any DFTTYP):}
\begin{shade}
convert4qmc -gamessAscii single det.out
\end{shade}
\begin{itemize}
\item{single det.out is the standard output generated by GAMESS.}
\end{itemize}
\item{\textit{(This option is not recommended. Use option below to avoid mistakes.)} For
multi-determinant calculations where the orbitals and configurations are read from different
files (for example when using orbitals from a MCSCF run and configurations from a
\item{\textit{(This option is not recommended. Use the following option to avoid mistakes.)} For
multideterminant calculations where the orbitals and configurations are read from different
files (e.g., when using orbitals from a MCSCF run and configurations from a
subsequent CI run):}
\begin{shade}
convert4qmc -gamessAscii orbitals multidet.out -ci cicoeff multidet.out
@ -660,7 +656,7 @@ convert4qmc -gamessAscii orbitals multidet.out -ci cicoeff multidet.out
orbitals. cicoeff multidet.out is the standard output from the calculation that calculates
the CI expansion.}
\end{itemize}
\item{For multi-determinant calculations where the orbitals and configurations are read from
\item{For multideterminant calculations where the orbitals and configurations are read from
the same file, using PRTMO=.T. in the GUESS input block:}
\begin{shade}
convert4qmc -gamessAscii multi det.out -ci multi det.out -readInitialGuess Norb
@ -673,17 +669,16 @@ convert4qmc -gamessAscii multi det.out -ci multi det.out -readInitialGuess Norb
Options:
\begin{itemize}
\item{\textbf{-gamessAscii file.out}: Standard output of GAMESS calculation. With the exception
of determinant configurations and coefficients in multi-determinant calculations,
everything else is read from this file including: atom coordinates, basis sets, single
particle orbitals, ECPs, number of electrons, multiplicity, etc.}
of determinant configurations and coefficients in multideterminant calculations,
everything else is read from this file including atom coordinates, basis sets, SPOs, ECPs, number of electrons, multiplicity, etc.}
\item{\textbf{-ci file.out}: In multi-determinant calculations, determinant configurations and
coefficients are read from this file. Notice that single particle orbitals are NOT read
from this file. Recognized CI packages are: ALDET, GUGA and ORMAS. Output
\item{\textbf{-ci file.out}: In multideterminant calculations, determinant configurations and
coefficients are read from this file. Notice that SPOs are NOT read
from this file. Recognized CI packages are ALDET, GUGA, and ORMAS. Output
produced with the GUGA package MUST have the option “NPRT=2” in the CIDRT
or DRT input blocks.}
\item{\textbf{-threshold cutoff}: Cutoff in multi-determinant expansion. Only configurations with
\item{\textbf{-threshold cutoff}: Cutoff in multideterminant expansion. Only configurations with
coefficients above this value are printed.}
\item{\textbf{-zeroCI}: Sets to zero the CI coefficients of all determinants, with the exception of the
@ -699,16 +694,16 @@ option in all CI calculations.}
\item{\textbf{-NaturalOrbitals Norb}: Read Norb “NATURAL ORBITALS” from GAMESS
output. The natural orbitals must exists in the output, otherwise the code aborts.}
\item{\textbf{-add3BodyJ}: Adds three-body Jastrow terms (e-e-I) between electron pairs (both
\item{\textbf{-add3BodyJ}: Adds 3-body Jastrow terms (e-e-I) between electron pairs (both
same spin and opposite spin terms) and all ion species in the system. The radial
function is initialized to zero and the default cutoff is 10.0 bohr. The converter will
add a one- and two-body Jastrow to the wavefunction block by default.}
function is initialized to zero, and the default cutoff is 10.0 bohr. The converter will
add a 1- and 2-body Jastrow to the wavefunction block by default.}
\end{itemize}
\subsection{Useful notes}
\begin{itemize}
\item{The type of single particle orbitals read by the converter depends on the type of
calculation and on the options used. By default, when neither -readInitialGuess or
\item{The type of SPOs read by the converter depends on the type of
calculation and on the options used. By default, when neither -readInitialGuess nor
-NaturalOrbitals are used, the following orbitals are read in each case (notice that
-readInitialGuess or -NaturalOrbitals are mutually exclusive):}
\begin{itemize}
@ -716,27 +711,27 @@ calculation and on the options used. By default, when neither -readInitialGuess
\item{MCSCF: “MCSCF OPTIMIZED ORBITALS”}
\item{GUGA, ALDET, ORMAS: Cannot read orbitals without -readInitialGuess or -NaturalOrbitals options.}
\end{itemize}
\item{The single particle orbitals and printed CI coefficients in MCSCF calculations are
\item{The SPOs and printed CI coefficients in MCSCF calculations are
not consistent in GAMESS. The printed CI coefficients correspond to the next-to-last
iteration, they are not recalculated with the final orbitals. So in order to get appropriate
iteration; they are not recalculated with the final orbitals. So to get appropriate
CI coefficients from MCSCF calculations, a subsequent CI (no SCF) calculation
is needed to produce consistent orbitals. In principle, it is possible to read the orbitals
from the MCSCF output and the CI coefficients and configurations from the
output of the following CI calculations. This could lead to problems in principle, since
GAMESS will rotate initial orbitals by default in order to obtain an initial guess consistent
output of the following CI calculations. This could lead to problems in principle since
GAMESS will rotate initial orbitals by default to obtain an initial guess consistent
with the symmetry of the molecule. This last step is done by default and can
change the orbitals reported in the MCSCF calculation before the CI is performed.
In order to avoid this problem, it is highly recommended to use option \#3 above to
read all the information from the output of the CI calculation, this requires the use
To avoid this problem, we highly recommend using the preceding option \#3 to
read all the information from the output of the CI calculation; this requires the use
of “PRTMO=.T.” in the GUESS input block. Since the orbitals are printed after any
symmetry rotation, the resulting output will always be consistent.}
\end{itemize}
\newpage
\section{Appendix C: Wave-function Optimization XML block}
\section{Appendix C: Wavefunction optimization XML block}\label{sec:lab_adv_mol_opt_appendix}
\begin{lstlisting}[style=QMCPXML,caption="Sample XML optimization block.",label=lst:lam_xml_opt]
\begin{lstlisting}[style=QMCPXML,caption=``Sample XML optimization block.",label=lst:lam_xml_opt]
<loop max="10">
<qmc method="linear" move="pbyp" checkpoint="-1" gpu="no">
<parameter name="blocks"> 10 </parameter>
@ -764,21 +759,21 @@ symmetry rotation, the resulting output will always be consistent.}
Options:
\begin{itemize}
\item{bigchange: (default 50.0) largest parameter change allowed}
\item{bigchange: (default 50.0) Largest parameter change allowed}
\item{usebuffer: (default no) Save useful information during VMC}
\item{nonlocalpp: (default no) Include non-local energy on 1-D min}
\item{nonlocalpp: (default no) Include nonlocal energy on 1-D min}
\item{MinMethod: (default quartic) Method to calculate magnitude of parameter change
quartic: fit quartic polynomial to 4 values of the cost function obtained using reweighting
quartic: fit quartic polynomial to four values of the cost function obtained using reweighting
along chosen direction linemin: direct line minimization using reweighting rescale:
no 1-D minimization. Uses Umrigars suggestions.}
\item{stepsize: (default 0.25) step size in either quartic or linemin methods.}
\item{alloweddifference: (default 1e-4) Allowed increased in energy}
\item{exp0: (default -16.0) Initial value for stabilizer (shift to diagonal of H) Actual value
\item{stepsize: (default 0.25) Step size in either quartic or linemin methods.}
\item{alloweddifference: (default 1e-4) Allowed increase in energy}
\item{exp0: (default -16.0) Initial value for stabilizer (shift to diagonal of H). Actual value
of stabilizer is 10 exp0}
\item{nstabilizers: (default 3) Number of stabilizers to try}
\item{stabilizaterScale: (default 2.0) Increase in value of exp0 between iterations.}
\item{max its: (default 1) number of inner loops with same sample}
\item{minwalkers: (default 0.3) minimum value allowed for the ratio of effective samples
\item{max its: (default 1) Number of inner loops with same sample}
\item{minwalkers: (default 0.3) Minimum value allowed for the ratio of effective samples
to actual number of walkers in a reweighting step. The optimization will stop if the
effective number of walkers in any reweighting calculation drops below this value. Last
set of acceptable parameters are kept.}
@ -791,15 +786,15 @@ Recommendations:
\item{Set samples to equal to (\#threads)*blocks.}
\item{Set steps to 1. Use substeps to control correlation between samples.}
\item{For cases where equilibration is slow, increase both substeps and warmupsteps.}
\item{For hard cases (e.g. simultaneous optimization of long MSD and 3-Body J), set exp0
\item{For hard cases (e.g., simultaneous optimization of long MSD and 3-Body J), set exp0
to 0 and do a single inner iteration (max its=1) per sample of configurations.}
\end{itemize}
\newpage
\section{Appendix D: VMC and DMC XML block}
\section{Appendix D: VMC and DMC XML block}\label{sec:lab_adv_mol_vmcdmc_appendix}
\begin{lstlisting}[style=QMCPXML,caption="Sample XML blocks for VMC and DMC calculations.",label=lst:lam_xml_vmc_dmc]
\begin{lstlisting}[style=QMCPXML,caption=``Sample XML blocks for VMC and DMC calculations.",label=lst:lam_xml_vmc_dmc]
<qmc method="vmc" move="pbyp" checkpoint="-1">
<parameter name="useDrift">yes</parameter>
<parameter name="warmupsteps">100</parameter>
@ -822,19 +817,19 @@ to 0 and do a single inner iteration (max its=1) per sample of configurations.}
General Options:
\begin{itemize}
\item{\textbf{move}: (default ”walker”) Type of electron move. Options: ”pbyp” and ”walker”.}
\item{\textbf{checkpoint}: (default ”-1”) (If > 0) Generate checkpoint files with given frequency.
\item{\textbf{move}: (default ``walker”) Type of electron move. Options: ``pbyp” and ``walker.”}
\item{\textbf{checkpoint}: (default ``-1”) (If > 0) Generate checkpoint files with given frequency.
The calculations can be restarted/continued with the produced checkpoint files.}
\item{\textbf{useDrift}: (default ”yes”) Defines the sampling mode. useDrift = ”yes” will
\item{\textbf{useDrift}: (default ``yes”) Defines the sampling mode. useDrift = ``yes” will
use Langevin acceleration to sample the VMC and DMC distributions, while
useDrift=”no” will use random displacements in a box.}
useDrift=``no” will use random displacements in a box.}
\item{\textbf{warmupSteps}: (default 0) Number of steps warmup steps at the beginning of the
calculation. No output is produced for these steps.}
\item{\textbf{blocks}: (default 1) Number of blocks (outer loop).}
\item{\textbf{steps}: (default 1) Number of steps per blocks (middle loop).}
\item{\textbf{sub steps}: (default 1) Number of substeps per step (inner loop). During sub steps,
\item{\textbf{sub steps}: (default 1) Number of substeps per step (inner loop). During substeps,
the local energy is not evaluated in VMC calculations, which leads to faster execution.
In VMC calculations, set sub steps to the average autocorrelation time of the desired
In VMC calculations, set substeps to the average autocorrelation time of the desired
quantity.}
\item{\textbf{time step}: (default 0.1) Electronic time step in bohr.}
\item{\textbf{samples}: (default 0) Number of walker configurations saved during the current
@ -845,18 +840,18 @@ number of walkers in the calculation will be equal to walkers*(\# nodes).}
Options unique to DMC:
\begin{itemize}
\item{\textbf{targetWalkers}: (default \#walkers from previous calculation, e.g. VMC.) Sets the
\item{\textbf{targetWalkers}: (default \#walkers from previous calculation, e.g., VMC). Sets the
target number of walkers. The actual population of walkers will fluctuate around this
value. The walkers will be distributed across all the nodes in the calculation. On a
given node, the walkers are split across all the threads in the system.}
\item{\textbf{nonlocalmoves}: (default ”no”) Set to ”yes” to turns on the use of Casula’s T-moves.}
\item{\textbf{nonlocalmoves}: (default ``no”) Set to ``yes” to turns on the use of Casula’s T-moves.}
\end{itemize}
\newpage
\section{Appendix E: Wave-function XML block}
\section{Appendix E: Wavefunction XML block}\label{sec:lab_adv_mol_wf_appendix}
\begin{lstlisting}[style=QMCPXML,caption="Basic framework for a single determinant determinantset XML block.",label=lst:lam_xml_determinantset]
\begin{lstlisting}[style=QMCPXML,caption=``Basic framework for a single-determinant determinantset XML block.",label=lst:lam_xml_determinantset]
<wavefunction name="psi0" target="e">
<determinantset type="MolecularOrbital" name="LCAOBSet"
source="ion0" transform="yes">
@ -897,25 +892,24 @@ these files (with very few exceptions like the cutoff of CI coefficients and the
functions) since changes can lead to unexpected results.
A QMCPACK wavefunction XML block is a combination of a determinantset, which
contains the anti-symmetric part of the wave-function, and one or more jastrow blocks.
The syntax of the anti-symmetric block depends on whether the wave-function is a single
determinant or a multi-determinant expansion. Listing~\ref{lst:lam_xml_determinantset}
contains the antisymmetric part of the wavefunction and one or more Jastrow blocks.
The syntax of the antisymmetric block depends on whether the wavefunction is a single
determinant or a multideterminant expansion. Listing \ref{lst:lam_xml_determinantset}
shows the general structure of the
single determinant case. The determinantset block is composed of a basisset block, which
defines the atomic orbital basis set, and a slaterdeterminant block, which defines the single
particle orbitals and occupation numbers of the Slater determinant.
% Listing~\ref{lst:lam_xml_basisset}
single-determinant case. The determinantset block is composed of a basisset block, which
defines the atomic orbital basis set, and a slaterdeterminant block, which defines the SPOs and occupation numbers of the Slater determinant.
% Listing \ref{lst:lam_xml_basisset}
%shows a section
%of a basisset block for a gold atom. The structure of this block is rigid and should not
%be modified.
Listing~\ref{lst:lam_xml_slaterdeterminant} shows a (piece of a) sample of a
Listing \ref{lst:lam_xml_slaterdeterminant} shows a (piece of a) sample of a
slaterdeterminant block. The
slaterdeterminant block consists of 2 determinant blocks, one for each electron spin. The
parameter “size” in the determinant block refers to the number of single particle orbitals
present while the “size” parameter in the coefficient block refers to the number of atomic
basis functions per single particle orbital.
slaterdeterminant block consists of two determinant blocks, one for each electron spin. The
parameter ``size” in the determinant block refers to the number of SPOs
present while the ``size” parameter in the coefficient block refers to the number of atomic
basis functions per SPO.
\begin{lstlisting}[style=QMCPXML,caption="Sample XML block for the single slater determinant case.",label=lst:lam_xml_slaterdeterminant]
\begin{lstlisting}[style=QMCPXML,caption=``Sample XML block for the single Slater determinant case.",label=lst:lam_xml_slaterdeterminant]
<slaterdeterminant>
<determinant id="updet" size="5">
<occupation mode="ground"/>
@ -935,17 +929,17 @@ basis functions per single particle orbital.
-0.00000000000000e+00 -3.20000000000000e-05 1.60000000000000e-05 1.60000000000000e-05
-0.00000000000000e+00 -0.00000000000000e+00 -0.00000000000000e+00 7.00000000000000e-06
\end{lstlisting}
Listing~\ref{lam_xml_multideterminant} shows the general structure of the multi-determinant case.
Listing \ref{lam_xml_multideterminant} shows the general structure of the multideterminant case.
Similar to the
single determinant case, the determinantset must contain a basisset block. This definition is
identical to the one described above. In this case, the definition of the single particle orbitals
must be done independently from the definition of the determinant configurations, the latter
is done in the sposet block while the former is done on the multideterminant block. Notice
that 2 sposet sets must be defined, one for each electron spin. The name of reach sposet set
single-determinant case, the determinantset must contain a basisset block. This definition is
identical to the one described previously. In this case, the definition of the SPOs
must be done independently from the definition of the determinant configurations; the latter
is done in the sposet block, while the former is done on the multideterminant block. Notice
that two sposet sets must be defined, one for each electron spin. The name of each sposet set
is required in the definition of the multideterminant block. The determinants are defined in
terms of occupation numbers based on these orbitals.
\begin{lstlisting}[style=QMCPXML,caption="Basic framework for a multi-determinant determinantset XML block.",label=lam_xml_multideterminant]
\begin{lstlisting}[style=QMCPXML,caption=``Basic framework for a multideterminant determinantset XML block.",label=lam_xml_multideterminant]
<wavefunction id="psi0" target="e">
<determinantset name="LCAOBSet" type="MolecularOrbital" transform="yes" source="ion0">
<basisset name="LCAOBSet">
@ -991,22 +985,22 @@ wavefunction optimization.}
\item{qchem coeff: (in csf ) Original coefficient of given configuration from GAMESS
calculation. This is used when applying a cutoff to the configurations read from the file.
The cutoff is applied on this parameter and not on the optimized coefficient.}
\item{nca and nab: number of core orbitals for up/down electrons. A core orbital is an
\item{nca and nab: Number of core orbitals for up/down electrons. A core orbital is an
orbital that is doubly occupied in all determinant configurations, not to be confused
with core electrons. These are not explicitly listed on the definition of configurations.}
\item{nea and neb: number of up/down active electrons (those being explicitly correlated).}
\item{nstates: number of correlated orbitals}
\item{size (in detlist ): contains the number of configurations in the list.}
\item{nea and neb: Number of up/down active electrons (those being explicitly correlated).}
\item{nstates: Number of correlated orbitals}.
\item{size (in detlist ): Contains the number of configurations in the list.}
\end{itemize}
The remaining part of the determinantset block is the definition of jastrow factor. Any
number of these can be defined. Figure~\ref{fig:lam_xml_jastrow} shows a sample jastrow
block including one-, two- and three-body terms. This is the standard block produced by
convert4qmc with the option -add3BodyJ (this particular example is for a water molecule).
The remaining part of the determinantset block is the definition of Jastrow factor. Any
number of these can be defined. Figure \ref{fig:lam_xml_jastrow} shows a sample Jastrow
block including 1-, 2- and 3-body terms. This is the standard block produced by
\texttt{convert4qmc} with the option -add3BodyJ (this particular example is for a water molecule).
Optimization of individual radial functions can be turned on/off using the “optimize”
parameter. It can be added to any coefficients block, even though it is currently not
present in the J1 and J2 blocks.
\begin{lstlisting}[style=QMCPXML,caption="Sample Jastrow XML block.",label=fig:lam_xml_jastrow]
\begin{lstlisting}[style=QMCPXML,caption=``Sample Jastrow XML block.",label=fig:lam_xml_jastrow]
<jastrow name="J2" type="Two-Body" function="Bspline" print="yes">
<correlation rcut="10" size="10" speciesA="u" speciesB="u">
<coefficients id="uu" type="Array">0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0</coefficients>
@ -1044,13 +1038,12 @@ present in the J1 and J2 blocks.
\end{lstlisting}
This training assumes basic familiarity with the UNIX operating system. In particular,
we use simple scripts written in “csh”. In addition, we assume that the student has obtained
all the necessary files and executables, and that the location of the training files are located
we use simple scripts written in “csh.” In addition, we assume you have obtained
all the necessary files and executables and that the training files are located
at \$\{TRAINING TOP\}.
The goal of the training not only to familiarize the student with the execution and
options in QMCPACK, but also to introduce him/her to important concepts in quantum
Monte Carlo calculations and many-body electronic structure calculations.
The goal of this training is not only to familiarize you with the execution and
options in QMCPACK but also to introduce you to important concepts in QMC calculations and many-body electronic structure calculations.

266
manual/lab_qmc_basics.tex
View File

@ -1,38 +1,38 @@
\chapter{Lab 2: QMC Basics}
\chapter{Lab 2: QMC basics}
\label{chap:lab_qmc_basics}
\section{Topics covered in this Lab}
This lab focuses on the basics of performing quality QMC calculations. As an example participants test an oxygen pseudopotential within DMC by calculating atomic and dimer properties, a common step prior to production runs. Topics covered include:
\section{Topics covered in this lab}
This lab focuses on the basics of performing quality QMC calculations. As an example, participants test an oxygen pseudopotential within DMC by calculating atomic and dimer properties, a common step prior to production runs. Topics covered include:
\begin{itemize}
\item{converting pseudopotentials into QMCPACK's FSATOM format}
\item{generating orbitals with Quantum ESPRESSO}
\item{converting orbitals into QMCPACK's ESHDF format with pw2qmcpack}
\item{optimizing Jastrow factors with QMCPACK}
\item{removing DMC timestep error via extrapolation}
\item{automating QMC workflows with Nexus}
\item{testing pseudopotentials for accuracy}
\item{Converting pseudopotentials into QMCPACK's FSATOM format}
\item{Generating orbitals with QE}
\item{Converting orbitals into QMCPACK's ESHDF format with pw2qmcpack}
\item{Optimizing Jastrow factors with QMCPACK}
\item{Removing DMC time step errors via extrapolation}
\item{Automating QMC workflows with Nexus}
\item{Testing pseudopotentials for accuracy}
\hide{
\item{(optional) running QMCPACK for a general system of interest}
\item{(optional) Running QMCPACK for a general system of interest}
}
\end{itemize}
\section{Lab outline}
\begin{enumerate}
\item{download and conversion of oxygen atom pseudopotential}
\item{DMC timestep study of the neutral oxygen atom}
\item{Download and conversion of oxygen atom pseudopotential}
\item{DMC time step study of the neutral oxygen atom}
\begin{enumerate}
\item{DFT orbital generation with Quantum ESPRESSO}
\item{orbital conversion with \ishell{pw2qmcpack.x}}
\item{optimization of Jastrow correlation factor with QMCPACK}
\item{DMC run with multiple timesteps}
\item{DFT orbital generation with QE}
\item{Orbital conversion with \ishell{pw2qmcpack.x}}
\item{Optimization of Jastrow correlation factor with QMCPACK}
\item{DMC run with multiple time steps}
\end{enumerate}
\item{DMC timestep study of the first ionization potential of oxygen}
\item{DMC time step study of the first ionization potential of oxygen}
\begin{enumerate}
\item{repetition of a-d above for ionized oxygen atom}
\item{Repetition of a-d above for ionized oxygen atom}
\end{enumerate}
\item{automated DMC calculations of the oxygen dimer binding curve}
\item{Automated DMC calculations of the oxygen dimer binding curve}
\end{enumerate}
@ -69,21 +69,21 @@ First enter the \ishell{oxygen_atom} directory:
cd labs/lab2_qmc_basics/oxygen_atom/
\end{shade}
\noindent
Throughout the rest of the lab, locations will be specified with respect to \ishell{labs/lab2_qmc_basics} (e.g. \ishell{oxygen_atom}).
Throughout the rest of the lab, locations are specified with respect to \ishell{labs/lab2_qmc_basics} (e.g., \ishell{oxygen_atom}).
We will use a potential from the Burkatzki-Filippi-Dolg pseudopotential database.
We use a potential from the Burkatzki-Filippi-Dolg pseudopotential database.
Although the full database is available in QMCPACK distribution (\ishell{trunk/pseudopotentials/BFD/}),
we use a BFD pseudopotential to illustrate the process of converting and testing an
external potential for use with QMCPACK. To obtain the pseudopotential, go to
\href{http://www.burkatzki.com/pseudos/index.2.html}{http://www.burkatzki.com/pseudos/index.2.html}
and click on the ``Select Pseudopotential'' button. Next click on oxygen in the
periodic table. Click on the empty circle next to ``V5Z'' (a large gaussian
basis set) and click on ``Next''. Select the Gamess format and click on
``Retrive Potential''. Helpful information about the pseudopotential will be
periodic table. Click on the empty circle next to ``V5Z'' (a large Gaussian
basis set) and click on ``Next.'' Select the Gamess format and click on
``Retrive Potential.'' Helpful information about the pseudopotential will be
displayed. The desired portion is at the bottom (the last 7 lines). Copy
this text into the editor of your choice (e.g. \ishell{emacs} or \ishell{vi})
this text into the editor of your choice (e.g., \ishell{emacs} or \ishell{vi})
and save it as \ishell{O.BFD.gamess}
(be sure to include a newline at the end of the file). To transform the
(be sure to include a new line at the end of the file). To transform the
pseudopotential into the FSATOM XML format used by QMCPACK, use the \ishell{ppconvert}
tool:
@ -101,9 +101,9 @@ jobrun_vesta ppconvert --gamess_pot O.BFD.gamess --s_ref "1s(2)2p(4)" \
\fi
\noindent
Observe the notation used to describe the reference valence configuration for this helium-core PP: \ishell{1s(2)2p(4)}. The \ishell{ppconvert} tool uses the following convention for the valence states: the first $s$ state is labeled \ishell{1s} (\ishell{1s}, \ishell{2s}, \ishell{3s}, \ldots), the first $p$ state is labeled \ishell{2p} (\ishell{2p}, \ishell{3p}, \ldots), the first $d$ state is labeled \ishell{3d} (\ishell{3d}, \ishell{4d}, \ldots). Copy the resulting xml file into the \ishell{oxygen_atom} directory.
Observe the notation used to describe the reference valence configuration for this helium-core PP: \ishell{1s(2)2p(4)}. The \ishell{ppconvert} tool uses the following convention for the valence states: the first $s$ state is labeled \ishell{1s} (\ishell{1s}, \ishell{2s}, \ishell{3s}, \ldots), the first $p$ state is labeled \ishell{2p} (\ishell{2p}, \ishell{3p}, \ldots), and the first $d$ state is labeled \ishell{3d} (\ishell{3d}, \ishell{4d}, \ldots). Copy the resulting xml file into the \ishell{oxygen_atom} directory.
Note: the command to convert the PP into QM Espresso's UPF format is similar (both formats are required):
Note: The command to convert the PP into QE's UPF format is similar (both formats are required):
\noindent
\ifws
@ -134,13 +134,13 @@ O-QMC GEN 2 1
The full QMCPACK pseudopotential is also included in \ishell{oxygen_atom/reference/O.BFD.*}.
\section{DFT with Quantum ESPRESSO to obtain the orbital part of the wavefunction}
\section{DFT with QE to obtain the orbital part of the wavefunction}
\label{sec:lqb_dft}
With the pseudopotential in hand, the next step toward a QMC calculation is to obtain the Fermionic part of the wavefunction, in this case a single Slater determinant constructed from DFT-LDA orbitals for a neutral oxygen atom. If you had trouble with the pseudopotential conversion step, pre-converted pseudopotential files are located in the \ishell{oxygen_atom/reference} directory.
With the pseudopotential in hand, the next step toward a QMC calculation is to obtain the Fermionic part of the wavefunction, in this case a single Slater determinant constructed from DFT-LDA orbitals for a neutral oxygen atom. If you had trouble with the pseudopotential conversion step, preconverted pseudopotential files are located in the \ishell{oxygen_atom/reference} directory.
Quantum ESPRESSO input for the DFT-LDA ground state of the neutral oxygen atom can be found in \ishell{O.q0.dft.in} and also listing~\ref{lst:O_q0_dft} below. Setting \ishell{wf_collect=.true.} instructs Quantum Espresso to write the orbitals to disk at the end of the run. Option \ishell{wf_collect=.true.} may be a potential problem in large simulations, it is recommended to avoid it and use the converter pw2qmcpack in parallel, see details in Sec.~\ref{sec:pw2qmcpack}. Note that the plane-wave energy cutoff has been set to a reasonable value of 300 Ry here (\ishell{ecutwfc=300}). This value depends on the pseudopotentials used, and in general should be selected by running DFT$\rightarrow$(orbital conversion)$\rightarrow$VMC with increasing energy cutoffs until the lowest VMC total energy and variance is reached.
QE input for the DFT-LDA ground state of the neutral oxygen atom can be found in \ishell{O.q0.dft.in} and also in listing \ref{lst:O_q0_dft}. Setting \ishell{wf_collect=.true.} instructs QE to write the orbitals to disk at the end of the run. Option \ishell{wf_collect=.true.} could be a potential problem in large simulations; therefore, we recommend avoiding it and using the converter pw2qmcpack in parallel (see details in Section~\ref{sec:pw2qmcpack}). Note that the plane-wave energy cutoff has been set to a reasonable value of 300 Ry here (\ishell{ecutwfc=300}). This value depends on the pseudopotentials used, and, in general, should be selected by running DFT$\rightarrow$(orbital conversion)$\rightarrow$VMC with increasing energy cutoffs until the lowest VMC total energy and variance is reached.
\begin{lstlisting}[style=ESPRESSO, language=espresso, caption={Quantum ESPRESSO input file for the neutral oxygen atom (\ishell{O.q0.dft.in})\label{lst:O_q0_dft}}]
\begin{lstlisting}[style=ESPRESSO, language=espresso, caption={QE input file for the neutral oxygen atom (\ishell{O.q0.dft.in})\label{lst:O_q0_dft}}.]
&CONTROL
calculation = 'scf'
restart_mode = 'from_scratch'
@ -192,7 +192,7 @@ CELL_PARAMETERS cubic
0.00000000 0.00000000 18.89726133
\end{lstlisting}
Run Quantum ESPRESSO by typing
Run QE by typing
\ifws
\begin{shade}
mpirun -np 4 pw.x -input O.q0.dft.in >&O.q0.dft.out&
@ -203,8 +203,7 @@ jobrun_vesta pw.x O.q0.dft.in
\end{shade}
\fi
The DFT run should take a few minutes to complete. If desired, you can track the progress of the DFT run by typing ``\ifws\ishell{tail -f O.q0.dft.out}
\else\ishell{tail -f O.q0.dft.output}\fi''. Once finished, you should check the LDA total energy in \ifws\ishell{O.q0.dft.out}\else\ishell{O.q0.dft.output}\fi by typing ``\ifws\ishell{grep '! ' O.q0.dft.out}\else\ishell{grep '! ' O.q0.dft.output}\fi''. The result should be close to
The DFT run should take a few minutes to complete. If desired, you can track the progress of the DFT run by typing ``\ifws\ishell{tail -f O.q0.dft.out}.\else\ishell{tail -f O.q0.dft.output}\fi'' Once finished, you should check the LDA total energy in \ifws\ishell{O.q0.dft.out}\: \else\ishell{O.q0.dft.output}\fi by typing ``\ifws\ishell{grep '! ' O.q0.dft.out}.\else\ishell{grep '! ' O.q0.dft.output}\fi'' The result should be close to
\begin{shade}
! total energy = -31.57553905 Ry
\end{shade}
@ -216,9 +215,9 @@ The DFT run should take a few minutes to complete. If desired, you can track th
The orbitals have been written in a format native to Quantum ESPRESSO in the \ishell{O.q0.save} directory. We will convert them into the ESHDF format expected by QMCPACK by using the \ishell{pw2qmcpack.x} tool. The input for \ishell{pw2qmcpack.x} can be found in the file \ishell{O.q0.p2q.in} and also in listing~\ref{lst:O_q0_p2q} below.
The orbitals have been written in a format native to QE in the \ishell{O.q0.save} directory. We will convert them into the ESHDF format expected by QMCPACK by using the \ishell{pw2qmcpack.x} tool. The input for \ishell{pw2qmcpack.x} can be found in the file \ishell{O.q0.p2q.in} and also in listing \ref{lst:O_q0_p2q}.
\begin{lstlisting}[caption={\ishell{pw2qmcpack.x} input file for orbital conversion (\ishell{O.q0.p2q.in})\label{lst:O_q0_p2q}}]
\begin{lstlisting}[caption={\ishell{pw2qmcpack.x} input file for orbital conversion (\ishell{O.q0.p2q.in})\label{lst:O_q0_p2q}}.]
&inputpp
prefix = 'O.q0'
outdir = './'
@ -241,10 +240,10 @@ Upon completion of the run, a new file should be present containing the orbitals
\section{Optimization with QMCPACK to obtain the correlated part of the wavefunction}\label{sec:optimization_walkthrough}
The wavefunction we have obtained to this point corresponds to a non-interacting Hamiltonian. Once the Coulomb pair potential is switched on between particles, it is known analytically that the exact wavefunction has cusps whenever two particles meet spatially and in general the electrons become correlated. This is represented in the wavefunction by introducing a Jastrow factor containing at least pair correlations
The wavefunction we have obtained to this point corresponds to a noninteracting Hamiltonian. Once the Coulomb pair potential is switched on between particles, it is known analytically that the exact wavefunction has cusps whenever two particles meet spatially and, in general, the electrons become correlated. This is represented in the wavefunction by introducing a Jastrow factor containing at least pair correlations:
\begin{align}
&\Psi_{Slater-Jastrow}=e^{-J}\Psi_{Slater} \\
&J = \sum_{\sigma\sigma'}\sum_{i<j}u^{\sigma\sigma'}_2(|r_i-r_j|) + \sum_\sigma\sum_{iI}u^{\sigma I}_1(|r_i-r_I|)
&J = \sum_{\sigma\sigma'}\sum_{i<j}u^{\sigma\sigma'}_2(|r_i-r_j|) + \sum_\sigma\sum_{iI}u^{\sigma I}_1(|r_i-r_I|)\:.
\end{align}
Here $\sigma$ is a spin variable while $r_i$ and $r_I$ represent electron and ion coordinates, respectively. The introduction of $J$ into the wavefunction is similar to F12 methods in quantum chemistry, though it has been present in essentially all QMC studies since the first applications the method (circa 1965).
@ -260,17 +259,17 @@ standard Slater-Jastrow form:
\Psi_T = e^{-(J_1+J_2)}D^\uparrow(\{\phi_u^\uparrow\}_{u=1}^{N^\uparrow})D^\downarrow(\{\phi_d^\downarrow\}_{d=1}^{N^\uparrow})
\end{align}
The orbitals forming the spin-restricted Slater determinants
($D^\uparrow/D^\downarrow$) are obtained from DFT or Hartree-Fock (\emph{e.g.} via Quantum ESPRESSO)
($D^\uparrow/D^\downarrow$) are obtained from DFT or Hartree-Fock (\emph{e.g.,} via QE)
and are fixed. The ground state of the (pseudo) oxygen atom is spin polarized
with $N^{\uparrow}=4$ and $N^{\downarrow}=2$.
The part of the wavefunction we will be optimizing is the Jastrow factor
($e^{-(J_1+J_2)}$), which in this case includes one- (electron-ion) and two-
($e^{-(J_1+J_2)}$), which in this case includes 1- (electron-ion) and 2-
(electron-electron) body correlation functions. The Jastrow factor is symmetric
under same-spin electron exchange and does not affect the DMC fixed node
approximation. Optimization of the Jastrow factor does, however, improve the
efficiency of the DMC calculation and reduces additional approximations due to
non-local pseudopotentials (locality approximation, T-moves). Note that a three-body
nonlocal pseudopotentials (locality approximation, T-moves). Note that a 3-body
term ($J_3$) is also available and is often necessary when using pseudopotentials
for transition metal species or when high accuracy is desired for molecules.
@ -284,19 +283,19 @@ for transition metal species or when high accuracy is desired for molecules.
}
\end{figure}
The explicit form of the one-body Jastrow factor we will be using is
The explicit form of the 1-body Jastrow factor we will be using is
\begin{align}\label{eq:J1}
J_1 = \sum_{e=1}^{N^\uparrow+N^\downarrow}U_1^{\uparrow/\downarrow}(|r_e-r_O|)
\end{align}
where $r_e$ refers to the electron positions and $r_O$ is
the position of the oxygen ion. The $U_1^{\uparrow/\downarrow}$ term is a
one-dimensional radial function represented with piecewise continuous cubic
1D radial function represented with piecewise continuous cubic
polynomials (B-splines). The adjustable parameters to be optimized are the
``knots'' of the B-splines which are simply the values of the $U_1$ function at
uniformly spaced grid points (See fig.~\ref{fig:u1_spline} for an example of a $U_1$
``knots'' of the B-splines, which are simply the values of the $U_1$ function at
uniformly spaced grid points (See Figure \ref{fig:u1_spline} for an example of a $U_1$
spline function with 8 knots).
The two-body Jastrow factor is spin resolved ($r^\uparrow/r^\downarrow$ are up/down electron positions):
The 2-body Jastrow factor is spin resolved ($r^\uparrow/r^\downarrow$ are up/down electron positions):
\begin{align}\label{eq:J2}
J_2 = \sum_{u<u'}U_2^{\uparrow\uparrow/\downarrow\downarrow}(|r_u^\uparrow-r_{u'}^\uparrow|) + \sum_{d<d'}U_2^{\uparrow\uparrow/\downarrow\downarrow}(|r_d^\downarrow-r_{d'}^\downarrow|) + \sum_{u,d} U_2^{\uparrow\downarrow}(|r_u^\uparrow-r_d^\downarrow|)
\end{align}
@ -311,14 +310,14 @@ Wavefunction optimization essentially relies on two inequalities regarding energ
E_T(P) &= \frac{\expvalh{\Psi_T(P)}{H}{\Psi_T(P)}}{\overlap{\Psi_T(P)}{\Psi_T(P)}} \ge E_0 \\
V_T(P) &= \frac{\expval{\Psi_T(P)}{\hat{H}^2}{\Psi_T(P)}}{\overlap{\Psi_T(P)}{\Psi_T(P)}} - \left(\frac{\expval{\Psi_T(P)}{H}{\Psi_T(P)}}{\overlap{\Psi_T(P)}{\Psi_T(P)}}\right)^2 \ge 0
\end{align}
Here $E_0$ is the ground state energy, $E_T(P)$ is the trial energy, $V_T(P)$ is the trial variance, and $P$ denotes the set of adjustable parameters in the trial wavefunction. Equality is reached only for the true ground state wavefunction and so the trial wavefunction can be improved by attempting to minimize a chosen cost function:
Here $E_0$ is the ground state energy, $E_T(P)$ is the trial energy, $V_T(P)$ is the trial variance, and $P$ denotes the set of adjustable parameters in the trial wavefunction. Equality is reached only for the true ground state wavefunction, and so the trial wavefunction can be improved by attempting to minimize a chosen cost function:
\begin{align}
C(P) = \alpha E_T(P) + (1-\alpha) V_T(P).
\end{align}
Iterative varational Monte Carlo methods have been developed to handle the non-linear optimization problem $\min\limits_P C(P)$. We will be using the linearized optimization method of Umrigar, \emph{et al.} (PRL \textbf{98} 110201 (2007)). Let us try this now with QMCPACK.
Iterative varational MC methods have been developed to handle the nonlinear optimization problem $\min\limits_P C(P)$. We will be using the linearized optimization method of Umrigar, \emph{et al.} (PRL \textbf{98} 110201 (2007)). Let us try this now with QMCPACK.
}
First, we need to update the template particle and wavefunction information in \ishell{O.q0.ptcl.xml} and \ishell{O.q0.wfs.xml}. We want to simulate the O atom in open boundary conditions (the default is periodic). To do this open \ishell{O.q0.ptcl.xml} with your favorite text editor (e.g. \ishell{emacs} or \ishell{vi}) and replace
First, we need to update the template particle and wavefunction information in \ishell{O.q0.ptcl.xml} and \ishell{O.q0.wfs.xml}. We want to simulate the O atom in open boundary conditions (the default is periodic). To do this, open \ishell{O.q0.ptcl.xml} with your favorite text editor (e.g., \ishell{emacs} or \ishell{vi}) and replace
\begin{lstlisting}[style=QMCPXML]
<parameter name="bconds">
p p p
@ -345,7 +344,7 @@ Next we will select Jastrow factors appropriate for an atom. In open boundary c
...
\end{lstlisting}
\noindent
These terms correspond to $u_2^{\uparrow\uparrow}/u_2^{\downarrow\downarrow}$, $u_2^{\uparrow\downarrow}$, and $u_1^{\uparrow O}/u_1^{\downarrow O}$, respectively. In each case, the correlation function ($u_*$) is represented by piecewise continuous cubic B-splines. Each correlation function has eight parameters which are just the values of $u$ on a uniformly spaced grid up to \ishell{rcut}. Initially the parameters (\ishell{coefficients}) are set to zero:
These terms correspond to $u_2^{\uparrow\uparrow}/u_2^{\downarrow\downarrow}$, $u_2^{\uparrow\downarrow}$, and $u_1^{\uparrow O}/u_1^{\downarrow O}$, respectively. In each case, the correlation function ($u_*$) is represented by piecewise continuous cubic B-splines. Each correlation function has eight parameters, which are just the values of $u$ on a uniformly spaced grid up to \ishell{rcut}. Initially the parameters (\ishell{coefficients}) are set to zero:
\begin{lstlisting}[style=QMCPXML]
<correlation speciesA="u" speciesB="u" size="8" rcut="10.0">
<coefficients id="uu" type="Array">
@ -387,20 +386,20 @@ The relevant portion of the input describing the linear optimization process is
</loop>
\end{lstlisting}
\noindent
An explanation of each input variable can be found below. The remaining variables control specialized internal details of the linear optimization algorithm. The meaning of these inputs is beyond the scope of this lab and reasonable results are often obtained keeping these values fixed.
An explanation of each input variable follows. The remaining variables control specialized internal details of the linear optimization algorithm. The meaning of these inputs is beyond the scope of this lab, and reasonable results are often obtained keeping these values fixed.
\begin{description}
\item[energy] Fraction of trial energy in the cost function.
\item[unreweightedvariance] Fraction of unreweighted trial variance in the cost function. Neglecting the weights can be more robust.
\item[reweightedvariance] Fraction of trial variance (including the full weights) in the cost function.
\item[timestep] Timestep of the VMC random walk, determines spatial distance moved by each electron during MC steps. Should be chosen such that the acceptance ratio of MC moves is around 50\% (30-70\% is often acceptable). Reasonable values are often between 0.2 and 0.6 $\textrm{Ha}^{-1}$.
\item[samples] Total number of MC samples collected for optimization, determines statistical error bar of cost function. Often efficient to start with a modest number of samples (50k) and then increase as needed. More samples may be required if the wavefunction contains a large number of variational parameters. MUST be be a multiple of the number of threads/cores \labsw{}{(use multiples of 512 on Vesta)}.
\item[timestep] Time step of the VMC random walk, determines spatial distance moved by each electron during MC steps. Should be chosen such that the acceptance ratio of MC moves is around 50\% (30--70\% is often acceptable). Reasonable values are often between 0.2 and 0.6 $\textrm{Ha}^{-1}$.
\item[samples] Total number of MC samples collected for optimization; determines statistical error bar of cost function. It is often efficient to start with a modest number of samples (50k) and then increase as needed. More samples may be required if the wavefunction contains a large number of variational parameters. MUST be be a multiple of the number of threads/cores\labsw{}{(use multiples of 512 on Vesta)}.
\item[warmupSteps] Number of MC steps discarded as a warmup or equilibration period of the random walk. If this is too small, it will bias the optimization procedure.
\item[blocks] Number of average energy values written to output files. Should be greater than 200 for meaningful statistical analysis of output data (\emph{e.g.} via \ishell{qmca}).
\item[subSteps] Number of MC steps in between energy evaluations. Each energy evaluation is expensive so taking a few steps to decorrelate between measurements can be more efficient. Will be less efficient with many substeps.
\item[nonlocalpp,useBuffer] If \ishell{nonlocalpp="no"}, then nonlocal part of the pseudopotential is not included when computing the cost function. If \ishell{useBuffer="yes"}, then temporary data is stored to speed up nonlocal pseudopotential evaluation at the expense of memory consumption.
\item[blocks] Number of average energy values written to output files. Should be greater than 200 for meaningful statistical analysis of output data ({e.g.,} via \ishell{qmca}).
\item[subSteps] Number of MC steps in between energy evaluations. Each energy evaluation is expensive, so taking a few steps to decorrelate between measurements can be more efficient. Will be less efficient with many substeps.
\item[nonlocalpp,useBuffer] If \ishell{nonlocalpp="no,"} then the nonlocal part of the pseudopotential is not included when computing the cost function. If \ishell{useBuffer="yes,"} then temporary data is stored to speed up nonlocal pseudopotential evaluation at the expense of memory consumption.
\item[loop max] Number of times to repeat the optimization. Using the resulting wavefunction from the previous optimization in the next one improves the results. Typical choices range between 8 and 16.
\end{description}
The cost function defines the quantity to be minimized during optimization. The three components of the cost function, energy, unreweighted variance, and reweighted variance should sum to one. Dedicating 100\% of the cost function to unreweighted variance is often a good choice. Another common choice is to try 90/10 or 80/20 mixtures of reweighted variance and energy. Using 100\% energy minimization is desirable for reducing DMC pseudopotential localization errors, but the optimization process is less stable and should only be attempted after performing several cycles of e.g. variance minimization first (the entire \ishell{loop} section can be duplicated with a different cost function each time).
The cost function defines the quantity to be minimized during optimization. The three components of the cost function, energy, unreweighted variance, and reweighted variance should sum to one. Dedicating 100\% of the cost function to unreweighted variance is often a good choice. Another common choice is to try 90/10 or 80/20 mixtures of reweighted variance and energy. Using 100\% energy minimization is desirable for reducing DMC pseudopotential localization errors, but the optimization process is less stable and should be attempted only after first performing several cycles of, for example, variance minimization (the entire \ishell{loop} section can be duplicated with a different cost function each time).
Replace \ishell{MAX}, \ishell{EVCOST}, \ishell{UVCOST}, \ishell{RVCOST}, \ishell{TS}, and \ishell{SAMPLES} in the \ishell{loop} with appropriate starting values in the \ishell{O.q0.opt.in.xml} input file. Perform the optimization run by typing
\ifws
@ -413,9 +412,9 @@ jobrun_vesta qmcpack O.q0.opt.in.xml
\end{shade}
\fi
\noindent
The run should only take a few minutes for reasonable values of loop \ishell{max} and \ishell{samples}.
The run should take only a few minutes for reasonable values of loop \ishell{max} and \ishell{samples}.
Log file output will appear in \labsw{\ishell{O.q0.opt.out}}{\ishell{O.q0.opt.output}}. The beginning of each linear optimization will be marked with text similar to
The log file output will appear in \labsw{\ishell{O.q0.opt.out}}{\ishell{O.q0.opt.output}}. The beginning of each linear optimization will be marked with text similar to
\begin{shade}
=========================================================
Start QMCFixedSampleLinearOptimize
@ -423,7 +422,7 @@ Log file output will appear in \labsw{\ishell{O.q0.opt.out}}{\ishell{O.q0.opt.ou
=========================================================
\end{shade}
\noindent
At the end of each optimization section the change in cost function, new values for the Jastrow parameters, and elapsed wallclock time are reported:
At the end of each optimization section the change in cost function, new values for the Jastrow parameters, and elapsed wall clock time are reported:
\begin{shade}
OldCost: 7.0598901869e-01 NewCost: 7.0592576381e-01 Delta Cost:-6.3254886314e-05
...
@ -457,7 +456,7 @@ eO_7 -4.1358075061e-03 1 1 ON 23
QMC Execution time = 2.8218972974e+01 secs
\end{shade}
\noindent
The cost function should decrease during each linear optimization (\ishell{Delta cost < 0}). Try ``\labsw{\ishell{grep OldCost *opt.out}}{\ishell{grep OldCost *opt.output}}''. You should see something like this:
The cost function should decrease during each linear optimization (\ishell{Delta cost < 0}). Try ``\labsw{\ishell{grep OldCost *opt.out.}}{\ishell{grep OldCost *opt.output}}'' You should see something like this:
\begin{shade}
OldCost: 1.2655186572e+00 NewCost: 7.2443875597e-01 Delta Cost:-5.4107990118e-01
OldCost: 7.2229830632e-01 NewCost: 6.9833678217e-01 Delta Cost:-2.3961524143e-02
@ -473,9 +472,9 @@ The cost function should decrease during each linear optimization (\ishell{Delta
OldCost: 7.0598901869e-01 NewCost: 7.0592576381e-01 Delta Cost:-6.3254886314e-05
\end{shade}
Blocked averages of energy data, including the kinetic energy and components of the potential energy, are written to \ishell{scalar.dat} files. The first is named ``\ishell{O.q0.opt.s000.scalar.dat}'', with a series number of zero (\ishell{s000}). In the end there will be \ishell{MAX} of them, one for each series.
Blocked averages of energy data, including the kinetic energy and components of the potential energy, are written to \ishell{scalar.dat} files. The first is named ``\ishell{O.q0.opt.s000.scalar.dat},'' with a series number of zero (\ishell{s000}). In the end there will be \ishell{MAX} of them, one for each series.
When the job has finished, use the \ishell{qmca} tool to assess the effectiveness of the optimization process. To look at just the total energy and the variance, type ``\ishell{qmca -q ev O.q0.opt*scalar*}''. This will print the energy, variance, and the variance/energy ratio in Hartree units:
When the job has finished, use the \ishell{qmca} tool to assess the effectiveness of the optimization process. To look at just the total energy and the variance, type ``\ishell{qmca -q ev O.q0.opt*scalar*}.'' This will print the energy, variance, and the variance/energy ratio in Hartree units:
\begin{shade}
LocalEnergy Variance ratio
O.q0.opt series 0 -15.739585 +/- 0.007656 0.887412 +/- 0.010728 0.0564
@ -498,7 +497,7 @@ Identify which optimization series is the ``best'' according to your cost functi
\textbf{\underline{Questions and Exercises}}
\end{flushleft}
\begin{enumerate}
\item{What is the acceptance ratio of your optimization runs? (use ``\ishell{qmca -q ar O.q0.opt*scalar*}'') Do you expect the Monte Carlo sampling to be efficient?}
\item{What is the acceptance ratio of your optimization runs? (use ``\ishell{qmca -q ar O.q0.opt*scalar*}'') Do you expect the MC sampling to be efficient?}
\item{How do you know when the optimization process has converged?}
% \item{Why is the mean and the error of the variance sometimes large? Consider using \newline``\ishell{qmca -t -q ev O.q0.opt*scalar*}'' to investigate.}
\item{(optional) Optimization is sometimes sensitive to initial guesses of the parameters. If you have time, try varying the initial parameters, including the cutoff radius (\ishell{rcut}) of the Jastrow factors (remember to change \ishell{id} in the \ishell{<project/>} element). Do you arrive at a similar set of final Jastrow parameters? What is the lowest variance you are able to achieve?}
@ -506,17 +505,17 @@ Identify which optimization series is the ``best'' according to your cost functi
\section{DMC timestep extrapolation I: neutral O atom}
The diffusion Monte Carlo (DMC) algorithm contains two biases in addition to the fixed node and pseudopotential approximations that are important to control: timestep and population control bias. In this section we will focus on estimating and removing timestep bias from DMC calculations. The essential fact to remember is that the bias vanishes as the timestep goes to zero while the needed computer time increases inversely with the timestep.
\section{DMC timestep extrapolation I: neutral oxygen atom}
The DMC algorithm contains two biases in addition to the fixed node and pseudopotential approximations that are important to control: time step and population control bias. In this section we focus on estimating and removing time step bias from DMC calculations. The essential fact to remember is that the bias vanishes as the time step goes to zero, while the needed computer time increases inversely with the time step.
% background on timestep error should be covered elsewhere in the manual
% perhaps replace this with a brief formula of error (order tau^2) on total energy
\hide{
The following subsection briefly discusses the origin of timestep and population control biases in DMC and how they can be minimized or extrapolated away. As before, the second subsection contains the lab walkthrough with QMCPACK. By the end of the section, we will have a solid DMC estimate of the ground state energy of oxygen.
The following subsection briefly discusses the origin of time step and population control biases in DMC and how they can be minimized or extrapolated away. As before, the second subsection contains the lab walkthrough with QMCPACK. By the end of the section, we will have a solid DMC estimate of the ground state energy of oxygen.
\subsubsection{Background on timestep and population control bias}\label{sec:opt_background}
\subsubsection{Background on time step and population control bias}\label{sec:opt_background}
DMC improves over the VMC algorithm by projecting toward the true many-body electronic ground state of the system. The projection operator is the (importance sampled) imaginary time propagator, which is also known as the thermodynamic density matrix:
\begin{align}
\hat{\rho} = e^{-t\hat{H}}
@ -525,22 +524,22 @@ The direct action of the projection operator on a trial wavefunction in position
\begin{align}
\expval{R}{e^{-t\hat{H}}}{\Psi_T} = \int dR' \rho(R,R';t)\Psi_T(R')
\end{align}
cannot be calculated in a straightforward fashion since the analytic form of $\rho(R,R';t)=\expval{R}{\rho}{R'}$ is unknown. In order to make the algorithm computationally tractable, the finite time projection operator is expanded as a product of short-time projection operators
cannot be calculated in a straightforward fashion since the analytic form of $\rho(R,R';t)=\expval{R}{\rho}{R'}$ is unknown. To make the algorithm computationally tractable, the finite-time projection operator is expanded as a product of short-time projection operators
\begin{align}
\expval{R}{e^{-t{H}}}{\Psi_T} &= \expval{R}{e^{-\tau\hat{H}}e^{-\tau\hat{H}}\cdots e^{-\tau\hat{H}}}{\Psi_T}\\
&=\int dR_1dR_2\cdots dR_M \rho(R,R_1;\tau)\rho(R_1,R_2;\tau)\cdots\rho(R_{M-1},R_M;\tau)\Psi_T(R_M)
\end{align}
The advantage here is that reasonable approximations of the short time propagators are known. Common approximations have the form
The advantage here is that reasonable approximations of the short-time propagators are known. Common approximations have the form
\begin{align}
\rho(R,R';\tau) = e^{D(R,R';\tau)}e^{B(R,R';\tau)} + \mathcal{O}(\tau^2)
\end{align}
where $D(R,R';\tau)$ and $B(R,R';\tau)$ represent drift and branching terms, respectively. DMC results are biased for any finite timestep ($\tau$). The bias can be eliminated by extrapolating to zero timestep. In practice this is done by performing a series of runs with decreasing timesteps and then fitting the results.
where $D(R,R';\tau)$ and $B(R,R';\tau)$ represent drift and branching terms, respectively. DMC results are biased for any finite time step ($\tau$). The bias can be eliminated by extrapolating to zero time step. In practice this is done by performing a series of runs with decreasing time steps and then fitting the results.
The drift term can be sampled with standard Monte Carlo methods, while the branching term is incorporated as a weight assigned to each random walker. Instead of accumulating the weight, it is more efficient to ``branch'' each walker according to the weight, resulting in some walkers being deleted and others copied multiple times. If left uncontrolled, the walker population $(P)$ may vanish or diverge. A stable algorithm is obtained by adjusting the branching weight to preserve the overall number of walkers on average. Population control also biases the results, but usually to a lesser extent than timestep error (the bias is proportional to $1/P$). A common rule of thumb is to use at least a couple thousand walkers. This bias should be checked occasionally by performing runs with varying numbers of walkers.
The drift term can be sampled with standard MC methods, while the branching term is incorporated as a weight assigned to each random walker. Instead of accumulating the weight, it is more efficient to ``branch'' each walker according to the weight, resulting in some walkers being deleted and others copied multiple times. If left uncontrolled, the walker population $(P)$ may vanish or diverge. A stable algorithm is obtained by adjusting the branching weight to preserve the overall number of walkers on average. Population control also biases the results, but usually to a lesser extent than time step error (the bias is proportional to $1/P$). A common rule of thumb is to use at least a couple thousand walkers. This bias should be checked occasionally by performing runs with varying numbers of walkers.
}
In the same directory you used to perform wavefunction optimization (\ishell{oxygen_atom}) you will find a sample DMC input file for the neutral oxygen atom named \ishell{O.q0.dmc.in.xml}. Open this file in a text editor and note the differences from the optimization case. Wavefunction information is no longer included from \ishell{pw2qmcpack}, but instead should come from the optimization run:
In the same directory you used to perform wavefunction optimization (\ishell{oxygen_atom}) you will find a sample DMC input file for the neutral oxygen atom named \ishell{O.q0.dmc.in.xml}. Open this file in a text editor and note the differences from the optimization case. Wavefunction information is no longer included from \ishell{pw2qmcpack} but instead should come from the optimization run:
\begin{lstlisting}[style=QMCPXML]
<!-- OPT_XML is from optimization, e.g. O.q0.opt.s008.opt.xml -->
<include href="OPT_XML"/>
@ -548,24 +547,24 @@ In the same directory you used to perform wavefunction optimization (\ishell{oxy
\noindent
Replace ``\ishell{OPT_XML}'' with the \ishell{opt.xml} file corresponding to the best Jastrow parameters you found in the last section (this is a file name similar to \ishell{O.q0.opt.s008.opt.xml}).
The QMC calculation section at the bottom is also different. The linear optimization blocks have been replaced with XML describing a VMC run followed by DMC. The input keywords are described below.
The QMC calculation section at the bottom is also different. The linear optimization blocks have been replaced with XML describing a VMC run followed by DMC. Descriptions of the input keywords follow.
\begin{description}
\item[timestep] Timestep of the VMC/DMC random walk. In VMC choose a timestep corresponding to an acceptance ratio of about 50\%. In DMC the acceptance ratio is often above 99\%.
\item[timestep] Time step of the VMC/DMC random walk. In VMC choose a time step corresponding to an acceptance ratio of about 50\%. In DMC the acceptance ratio is often above 99\%.
\item[warmupSteps] Number of MC steps discarded as a warmup or equilibration period of the random walk.
\item[steps] Number of MC steps per block. Physical quantities, such as the total energy, are averaged over walkers and steps.
\item[blocks] Number of blocks. This is also the number of average energy values written to output files. Should be greater than 200 for meaningful statistical analysis of output data (\emph{e.g.} via \ishell{qmca}). The total number of MC steps each walker takes is \ishell{blocks}$\times$\ishell{steps}.
\item[blocks] Number of blocks. This is also the number of average energy values written to output files. The number should be greater than 200 for meaningful statistical analysis of output data (e.g., via \ishell{qmca}). The total number of MC steps each walker takes is \ishell{blocks}$\times$\ishell{steps}.
\item[samples] VMC only. This is the number of walkers used in subsequent DMC runs. Each DMC walker is initialized with electron positions sampled from the VMC random walk.
\item[nonlocalmoves] DMC only. If yes/no, use the locality approximation/T-moves for non-local pseudopotentials. T-moves generally improve the stability of the algorithm and restore the variational principle for small systems (T-moves version 1).
\item[nonlocalmoves] DMC only. If yes/no, use the locality approximation/T-moves for nonlocal pseudopotentials. T-moves generally improve the stability of the algorithm and restore the variational principle for small systems (T-moves version 1).
\end{description}
The purpose of the VMC run is to provide initial electron positions for each DMC walker. Setting $texttt{walkers}]=1$ in the VMC block ensures there will be only one VMC walker per execution thread. There will be a total of \labsw{4}{512} VMC walkers in this case (see \ishell{O.q0.dmc.qsub.in}). We want the electron positions used to initialize the DMC walkers to be decorrelated from one another. A VMC walker will often decorrelate from its current position after propagating for a few Ha$^{-1}$ in imaginary time (in general this is system dependent). This leads to a rough rule of thumb for choosing \ishell{blocks} and \ishell{steps} for the VMC run (\labsw{$\texttt{VWALKERS}=4$}{$\texttt{VWALKERS}=512$} here):
The purpose of the VMC run is to provide initial electron positions for each DMC walker. Setting $texttt{walkers}]=1$ in the VMC block ensures there will be only one VMC walker per execution thread. There will be a total of \labsw{4}{512} VMC walkers in this case (see \ishell{O.q0.dmc.qsub.in}). We want the electron positions used to initialize the DMC walkers to be decorrelated from one another. A VMC walker will often decorrelate from its current position after propagating for a few Ha$^{-1}$ in imaginary time (in general, this is system dependent). This leads to a rough rule of thumb for choosing \ishell{blocks} and \ishell{steps} for the VMC run (\labsw{$\texttt{VWALKERS}=4$}{$\texttt{VWALKERS}=512$} here):
\begin{align}
\texttt{VBLOCKS}\times\texttt{VSTEPS} \ge \frac{\texttt{DWALKERS}}{\texttt{VWALKERS}} \frac{5~\textrm{Ha}^{-1}}{\texttt{VTIMESTEP}}
\end{align}
Fill in the VMC XML block with appropriate values for these parameters. There should be more than one DMC walker per thread and enough walkers in total to avoid population control bias. The general rule of thumb is to have more than $\sim 2000$ walkers, although the dependence of the total energy on population size should be explicitly checked from time to time.
Fill in the VMC XML block with appropriate values for these parameters. There should be more than one DMC walker per thread and enough walkers in total to avoid population control bias. The general rule of thumb is to have more than $\sim 2,000$ walkers, although the dependence of the total energy on population size should be explicitly checked from time to time.
To study timestep bias, we will perform a sequence of DMC runs over a range of timesteps ($0.1$ Ha$^{-1}$ is too large and timesteps below $0.002$ Ha$^{-1}$ are probably too small). A common approach is to select a fairly large timestep to begin with and then decrease the timestep by a factor of two in each subsequent DMC run. The total amount of imaginary time the walker population propagates should be the same for each run. A simple way to accomplish this is to choose input parameters in the following way
To study time step bias, we will perform a sequence of DMC runs over a range of time steps ($0.1$ Ha$^{-1}$ is too large, and time steps below $0.002$ Ha$^{-1}$ are probably too small). A common approach is to select a fairly large time step to begin with and then decrease the time step by a factor of two in each subsequent DMC run. The total amount of imaginary time the walker population propagates should be the same for each run. A simple way to accomplish this is to choose input parameters in the following way
\begin{align}\label{eq:timestep_iter}
\texttt{timestep}_{n} &= \texttt{timestep}_{n-1}/2\nonumber\\
\texttt{warmupSteps}_{n} &= \texttt{warmupSteps}_{n-1}\times 2\nonumber\\
@ -574,7 +573,7 @@ To study timestep bias, we will perform a sequence of DMC runs over a range of t
\end{align}
Each DMC run will require about twice as much computer time as the one preceding it. Note that the number of blocks is kept fixed for uniform statistical analysis. $\texttt{blocks}\times\texttt{steps}\times\texttt{timestep}\sim 60~\mathrm{Ha}^{-1}$ is sufficient for this system.
Choose an initial DMC timestep and create a sequence of $N$ timesteps according to~\ref{eq:timestep_iter}. Make $N$ copies of the DMC XML block in the input file
Choose an initial DMC time step and create a sequence of $N$ time steps according to \ref{eq:timestep_iter}. Make $N$ copies of the DMC XML block in the input file.
\begin{lstlisting}[style=QMCPXML]
<qmc method="dmc" move="pbyp">
<parameter name="warmupSteps" > DWARMUP </parameter>
@ -585,7 +584,7 @@ Choose an initial DMC timestep and create a sequence of $N$ timesteps according
</qmc>
\end{lstlisting}
\noindent
Fill in \ishell{DWARMUP}, \ishell{DBLOCKS}, \ishell{DSTEPS}, and \ishell{DTIMESTEP} for each DMC run according to~\ref{eq:timestep_iter}. Start the DMC timestep extrapolation run by typing:
Fill in \ishell{DWARMUP}, \ishell{DBLOCKS}, \ishell{DSTEPS}, and \ishell{DTIMESTEP} for each DMC run according to \ref{eq:timestep_iter}. Start the DMC time step extrapolation run by typing:
\ifws
\begin{shade}
mpirun -np 4 qmcpack O.q0.dmc.in.xml >&O.q0.dmc.out&
@ -598,9 +597,9 @@ jobrun_vesta qmcpack O.q0.dmc.in.xml
\noindent
The run should take only a few minutes to complete.
QMCPACK will create files prefixed with \ishell{O.q0.dmc}. The log file is \labsw{\ishell{O.q0.dmc.out}}{\ishell{O.q0.dmc.output}}. As before, block averaged data is written to \ishell{scalar.dat} files. In addition, DMC runs produce \ishell{dmc.dat} files which contain energy data averaged only over the walker population (one line per DMC step). The \ishell{dmc.dat} files also provide a record of the walker population at each step.
QMCPACK will create files prefixed with \ishell{O.q0.dmc}. The log file is \labsw{\ishell{O.q0.dmc.out}}{\ishell{O.q0.dmc.output}}. As before, block-averaged data is written to \ishell{scalar.dat} files. In addition, DMC runs produce \ishell{dmc.dat} files, which contain energy data averaged only over the walker population (one line per DMC step). The \ishell{dmc.dat} files also provide a record of the walker population at each step.
Use the \ishell{PlotTstepConv.pl} to obtain a linear fit to the timestep data (type ``\ishell{PlotTstepConv.pl O.q0.dmc.in.xml 40}''). You should see a plot similar to fig.~\ref{fig:timestep_conv}. The tail end of the text output displays the parameters for the linear fit. The ``\ishell{a}'' parameter is the total energy extrapolated to zero timestep in Hartree units.
Use the \ishell{PlotTstepConv.pl} to obtain a linear fit to the time step data (type ``\ishell{PlotTstepConv.pl O.q0.dmc.in.xml 40}''). You should see a plot similar to Figure \ref{fig:timestep_conv}. The tail end of the text output displays the parameters for the linear fit. The ``\ishell{a}'' parameter is the total energy extrapolated to zero time step in Hartree units.
\begin{shade}
...
@ -631,22 +630,22 @@ b = -0.0457479 +/- 0.0422 (92.24%)
\end{flushleft}
\begin{enumerate}
\item{What is the $\tau\rightarrow 0$ extrapolated value for the total energy?}
\item{What is the maximum timestep you should use if you want to calculate the total energy to an accuracy of $0.05$ eV? For convenience, $1~\textrm{Ha}=27.2113846~\textrm{eV}$.}
\item{What is the acceptance ratio for this (bias$<0.05$ eV) run? Does it follow the rule of thumb for sensible DMC (acceptance ratio $>99$\%) ?}
\item{What is the maximum time step you should use if you want to calculate the total energy to an accuracy of $0.05$ eV? For convenience, $1~\textrm{Ha}=27.2113846~\textrm{eV}$.}
\item{What is the acceptance ratio for this (bias $<0.05$ eV) run? Does it follow the rule of thumb for sensible DMC (acceptance ratio $>99$\%) ?}
\item{Check the fluctuations in the walker population (\ishell{qmca -t -q nw O.q0.dmc*dmc.dat --noac}). Does the population seem to be stable?}
\item{(Optional) Study population control bias for the oxygen atom. Select a few population sizes \labsw{}{(use multiples of 512 to fit cleanly on a single Vesta partition)}. Copy \ishell{O.q0.dmc.in.xml} to a new file and remove all but one DMC run (select a single timestep). Make one copy of the new file for each population, set ``\ishell{samples}'', and choose a unique \ishell{id} in \ishell{<project/>}. \labsw{}{Run one job at a time to avoid crowding the lab allocation.} Use \ishell{qmca} to study the dependence of the DMC total energy on the walker population. How large is the bias compared to timestep error? What bias is incurred by following the ``rule of thumb'' of a couple thousand walkers? Will population control bias generally be an issue for production runs on modern parallel machines?}
\item{(Optional) Study population control bias for the oxygen atom. Select a few population sizes. \labsw{}{(use multiples of 512 to fit cleanly on a single Vesta partition)} Copy \ishell{O.q0.dmc.in.xml} to a new file and remove all but one DMC run (select a single time step). Make one copy of the new file for each population, set ``\ishell{samples},'' and choose a unique \ishell{id} in \ishell{<project/>}. \labsw{}{Run one job at a time to avoid crowding the lab allocation.} Use \ishell{qmca} to study the dependence of the DMC total energy on the walker population. How large is the bias compared with time step error? What bias is incurred by following the ``rule of thumb'' of a couple thousand walkers? Will population control bias generally be an issue for production runs on modern parallel machines?}
\end{enumerate}
\section{DMC timestep extrapolation II: O atom ionization potential}
In this section, we will repeat the calculations of the prior two sections (optimization, timestep extrapolation) for the $+1$ charge state of the oxygen atom. Comparing the resulting 1st ionization potential (IP) with experimental data will complete our first test of the BFD oxygen pseudopotential. In actual practice, higher IP's could also be tested prior to performing production runs.
\section{DMC time step extrapolation II: oxygen atom ionization potential}
In this section, we will repeat the calculations of the previous two sections (optimization, time step extrapolation) for the $+1$ charge state of the oxygen atom. Comparing the resulting first ionization potential (IP) with experimental data will complete our first test of the BFD oxygen pseudopotential. In actual practice, higher IPs could also be tested before performing production runs.
Obtaining the timestep extrapolated DMC total energy for ionized oxygen should take much less (human) time than for the neutral case. For convenience, the necessary steps are briefly summarized below.
Obtaining the time step extrapolated DMC total energy for ionized oxygen should take much less (human) time than for the neutral case. For convenience, the necessary steps are summarized as follows.
\begin{enumerate}
\item{Obtain DFT orbitals with Quantum ESPRESSO}
\item{Obtain DFT orbitals with QE.}
\begin{enumerate}
\item{Copy the DFT input (\ishell{O.q0.dft.in}) to \ishell{O.q1.dft.in}}
\item{Edit \ishell{O.q1.dft.in} to match the +1 charge state of the oxygen atom}
\item{Edit \ishell{O.q1.dft.in} to match the +1 charge state of the oxygen atom.}
\begin{lstlisting}[style=espresso]
...
prefix = 'O.q1'
@ -659,10 +658,10 @@ Obtaining the timestep extrapolated DMC total energy for ionized oxygen should t
\else\ishell{jobrun_vesta pw.x O.q1.dft.in}\fi}
\end{enumerate}
\item{Convert the orbitals to ESHDF format}
\item{Convert the orbitals to ESHDF format.}
\begin{enumerate}
\item{Copy the pw2qmcpack input (\ishell{O.q0.p2q.in}) to \ishell{O.q1.p2q.in}}
\item{Edit \ishell{O.q1.p2q.in} to match the file prefix used in DFT}
\item{Edit \ishell{O.q1.p2q.in} to match the file prefix used in DFT.}
\begin{verbatim}
...
prefix = 'O.q1'
@ -671,10 +670,10 @@ Obtaining the timestep extrapolated DMC total energy for ionized oxygen should t
\item{Perform the orbital conversion run: \ifws\ishell{mpirun -np 1 pw2qmcpack.x<O.q1.p2q.in>&O.q1.p2q.out&}\else\ishell{jobrun_vesta pw2qmcpack.x O.q1.p2q.in}\fi}
\end{enumerate}
\item{Optimize the Jastrow factor with QMCPACK}
\item{Optimize the Jastrow factor with QMCPACK.}
\begin{enumerate}
\item{Copy the optimization input (\ishell{O.q0.opt.in.xml}) to \ishell{O.q1.opt.in.xml}}
\item{Edit \ishell{O.q1.opt.in.xml} to match the file prefix used in DFT}
\item{Edit \ishell{O.q1.opt.in.xml} to match the file prefix used in DFT.}
\begin{lstlisting}[style=QMCPXML]
...
<project id="O.q1.opt" series="0">
@ -684,13 +683,13 @@ Obtaining the timestep extrapolated DMC total energy for ionized oxygen should t
<include href="O.q1.wfs.xml"/>
...
\end{lstlisting}
\item{Edit the particle XML file (O.q1.ptcl.xml) to have open boundary conditions}
\item{Edit the particle XML file (\ishell{O.q1.ptcl.xml}) to have open boundary conditions.}
\begin{lstlisting}[style=QMCPXML]
<parameter name="bconds">
n n n
</parameter>
\end{lstlisting}
\item{Add cutoffs to the Jastrow factors in the wavefunction XML file (O.q1.wfs.xml)}
\item{Add cutoffs to the Jastrow factors in the wavefunction XML file (\ishell{O.q1.wfs.xml}).}
\begin{lstlisting}[style=QMCPXML]
...
<correlation speciesA="u" speciesB="u" size="8" rcut="10.0">
@ -704,10 +703,10 @@ Obtaining the timestep extrapolated DMC total energy for ionized oxygen should t
\item{Identify the optimal set of parameters with \ishell{qmca} (\ishell{[your opt.xml]}).}
\end{enumerate}
\item{DMC timestep study with QMCPACK}
\item{DMC time step study with QMCPACK}
\begin{enumerate}
\item{Copy the DMC input (\ishell{O.q0.dmc.in.xml}) to \ishell{O.q1.dmc.in.xml}}
\item{Edit \ishell{O.q1.dmc.in.xml} to use the DFT prefix and the optimal Jastrow}
\item{Edit \ishell{O.q1.dmc.in.xml} to use the DFT prefix and the optimal Jastrow.}
\begin{lstlisting}[style=QMCPXML]
...
<project id="O.q1.dmc" series="0">
@ -718,23 +717,23 @@ Obtaining the timestep extrapolated DMC total energy for ionized oxygen should t
...
\end{lstlisting}
\item{Perform the DMC run: \ifws\ishell{mpirun -np 4 qmcpack O.q1.dmc.in.xml >&O.q1.dmc.out&}\else\ishell{jobrun_vesta qmcpack O.q1.dmc.in.xml}\fi}
\item{Obtain the DMC total energy extrapolated to zero timestep with \ishell{PlotTstepConv.pl}.}
\item{Obtain the DMC total energy extrapolated to zero time step with \ishell{PlotTstepConv.pl}.}
\end{enumerate}
\end{enumerate}
The process listed above, which excludes additional steps for orbital generation and conversion, can become tedious to perform by hand in production settings where many calculations are often required. For this reason automation tools are introduced for calculations involving the oxygen dimer in section~\ref{sec:dimer_automation} of the lab.
The aforementioned process, which excludes additional steps for orbital generation and conversion, can become tedious to perform by hand in production settings where many calculations are often required. For this reason, automation tools are introduced for calculations involving the oxygen dimer in Section \ref{sec:dimer_automation} of the lab.
\vspace{1cm}
\begin{flushleft}
\textbf{\underline{Questions and Exercises}}
\end{flushleft}
\begin{enumerate}
\item{What is the $\tau\rightarrow 0$ extrapolated DMC value for the 1st ionization potential of oxygen?}
\item{How does the extrapolated value compare to the experimental IP? Go to\newline \href{http://physics.nist.gov/PhysRefData/ASD/ionEnergy.html}{http://physics.nist.gov/PhysRefData/ASD/ionEnergy.html} and enter ``\ishell{O I}'' in the box labeled ``\ishell{Spectra}'' and click on the ``\ishell{Retrieve Data}'' button.
\item{What is the $\tau\rightarrow 0$ extrapolated DMC value for the first ionization potential of oxygen?}
\item{How does the extrapolated value compare with the experimental IP? Go to\newline \href{http://physics.nist.gov/PhysRefData/ASD/ionEnergy.html}{http://physics.nist.gov/PhysRefData/ASD/ionEnergy.html} and enter ``\ishell{O I}'' in the box labeled ``\ishell{Spectra}'' and click on the ``\ishell{Retrieve Data}'' button.
%For comparison the LDA value is $12.25$ eV.
}
\item{What can we conclude about the accuracy of the pseudopotential? What factors complicate this assessment?}
\item{Explore the sensitivity of the IP to the choice of timestep. Type ``\ishell{./ip_conv.py}'' to view three timestep extrapolation plots: two for the $q=0,1$ total energies and one for the IP. Is the IP more, less, or similarly sensitive to timestep than the total energy?}
\item{What is the maximum timestep you should use if you want to calculate the ionization potential to an accuracy of $0.05$ eV? What factor of cpu time is saved by assessing timestep convergence on the IP (a total energy difference) vs. a single total energy?}
\item{Explore the sensitivity of the IP to the choice of time step. Type ``\ishell{./ip_conv.py}'' to view three time step extrapolation plots: two for the $q=0,$ one for total energies, and one for the IP. Is the IP more, less, or similarly sensitive to time step than the total energy?}
\item{What is the maximum time step you should use if you want to calculate the ionization potential to an accuracy of $0.05$ eV? What factor of CPU time is saved by assessing time step convergence on the IP (a total energy difference) vs. a single total energy?}
\item{Are the acceptance ratio and population fluctuations reasonable for the $q=1$ calculations?}
\end{enumerate}
@ -744,16 +743,16 @@ The process listed above, which excludes additional steps for orbital generation
\section{DMC workflow automation with Nexus}
Production QMC projects are often composed of many similar workflows. The simplest of these is a single DMC calculation involving four different compute jobs:
\begin{enumerate}
\item{Orbital generation via Quantum ESPRESSO or GAMESS.}
\item{Orbital generation via QE or GAMESS.}
\item{Conversion of orbital data via \ishell{pw2qmcpack.x} or \ishell{convert4qmc}.}
\item{Optimization of Jastrow factors via QMCPACK.}
\item{DMC calculation via QMCPACK.}
\end{enumerate}
Simulation workflows quickly become more complex with increasing costs in terms of human time for the researcher. Automation tools can decrease both human time and error if used well.
The set of automation tools we will be using is known as Nexus \cite{Krogel2016nexus}, which is distributed with QMCPACK. Nexus is capable of generating input files, submitting and monitoring compute jobs, passing data between simulations (such as relaxed structures, orbital files, optimized Jastrow parameters, etc.), and data analysis. The user interface to Nexus is through a set of functions defined in the Python programming language. User scripts that execute simple workflows resemble input files and do not require programming experience. More complex workflows require only basic programming constructs (\emph{e.g.} for loops and if statements). Nexus input files/scripts should be easier to navigate than QMCPACK input files and more efficient than submitting all the jobs by hand.
The set of automation tools we will be using is known as Nexus \cite{Krogel2016nexus}, which is distributed with QMCPACK. Nexus is capable of generating input files, submitting and monitoring compute jobs, passing data between simulations (relaxed structures, orbital files, optimized Jastrow parameters, etc.), and data analysis. The user interface to Nexus is through a set of functions defined in the Python programming language. User scripts that execute simple workflows resemble input files and do not require programming experience. More complex workflows require only basic programming constructs (e.g. for loops and if statements). Nexus input files/scripts should be easier to navigate than QMCPACK input files and more efficient than submitting all the jobs by hand.
Nexus is driven by simple user-defined scripts that resemble keyword-driven input files. An example Nexus input file that performs a single VMC calculation (with pre-generated orbitals) is shown below. Take a moment to read it over and especially note the comments (prefixed with ``\ishell{\#}'') explaining most of the contents. If the input syntax is unclear you may want to consult portions of appendix~\ref{app:python_basics}, which gives a condensed summary of Python constructs. An additional example and details about the inner workings of Nexus can be found in the reference publication \cite{Krogel2016nexus}.
Nexus is driven by simple user-defined scripts that resemble keyword-driven input files. An example Nexus input file that performs a single VMC calculation (with pregenerated orbitals) follows. Take a moment to read it over and especially note the comments (prefixed with ``\ishell{\#}'') explaining most of the contents. If the input syntax is unclear you may want to consult portions of Appendix \ref{app:python_basics}, which gives a condensed summary of Python constructs. An additional example and details about the inner workings of Nexus can be found in the reference publication \cite{Krogel2016nexus}.
%For more information about the functionality and effective use of Nexus, consult \ishell{docs/Nexus.pdf}.
@ -966,7 +965,8 @@ The main difference is that a full workflow of runs (DFT orbital generation, orb
%\end{lstlisting}
As in the example in the last section, the oxygen dimer is generated with the \ishell{generate_physical_system} function:
As in the example in the last section, the oxygen dimer is generated with the \ishell{generate_physical_
system} function:
\begin{lstlisting}[style=Python]
dimer = generate_physical_system(
type = 'dimer',
@ -981,7 +981,7 @@ dimer = generate_physical_system(
\noindent
Similar syntax can be used to generate crystal structures or to specify systems with arbitrary atomic configurations and simulation cells. Notice that a ``\ishell{scale}'' variable has been introduced to stretch or compress the dimer.
Next, objects representing a Quantum ESPRESSO (PWSCF) run and subsequent orbital conversion step are constructed with respective \ishell{generate_*} functions:
Next, objects representing a QE (PWSCF) run and subsequent orbital conversion step are constructed with respective \ishell{generate_*} functions:
\begin{lstlisting}[style=Python]
dft = generate_pwscf(
identifier = 'dft',
@ -999,7 +999,7 @@ p2q = generate_pw2qmcpack(
)
sims.append(p2q)
\end{lstlisting}
Note the \ishell{dependencies} keyword. This keyword is used to construct workflows out of otherwise separate runs. In this case, the dependency indicates that the orbital conversion run must wait for the DFT to finish prior to starting.
Note the \ishell{dependencies} keyword. This keyword is used to construct workflows out of otherwise separate runs. In this case, the dependency indicates that the orbital conversion run must wait for the DFT to finish before starting.
Objects representing QMCPACK simulations are then constructed with the \ishell{generate_qmcpack} function:
\begin{lstlisting}[style=Python]
@ -1067,14 +1067,14 @@ qmc = generate_qmcpack(
sims.append(qmc)
\end{lstlisting}
\noindent
Shared details such as the run directory, job, pseudopotentials, and orbital file have been omitted (\ishell{...}). The ``\ishell{opt}'' run will optimize a 1-body B-spline Jastrow with 8 knots having a cutoff of 5.0 Bohr and a B-spline Jastrow (for up-up and up-down correlations) with 8 knots and cutoffs of 10.0 Bohr. The Jastrow list for the DMC run is empty and the usage of \ishell{dependencies} above indicates that the DMC run depends on the optimization run for the Jastrow factor. Nexus will submit the ``\ishell{opt}'' run first and upon completion it will scan the output, select the optimal set of parameters, pass the Jastrow information to the ``\ishell{qmc}'' run and then submit the DMC job. Independent job workflows are submitted in parallel when permitted \labsw{}{(we have explicitly limited this for this lab by setting \ishell{queue_size=2} for Vesta)}. No input files are written or job submissions made until the ``\ishell{run_project}'' function is reached:
Shared details such as the run directory, job, pseudopotentials, and orbital file have been omitted (\ishell{...}). The ``\ishell{opt}'' run will optimize a 1-body B-spline Jastrow with 8 knots having a cutoff of 5.0 Bohr and a B-spline Jastrow (for up-up and up-down correlations) with 8 knots and cutoffs of 10.0 Bohr. The Jastrow list for the DMC run is empty, and the previous use of \ishell{dependencies} indicates that the DMC run depends on the optimization run for the Jastrow factor. Nexus will submit the ``\ishell{opt}'' run first, and upon completion it will scan the output, select the optimal set of parameters, pass the Jastrow information to the ``\ishell{qmc}'' run, and then submit the DMC job. Independent job workflows are submitted in parallel when permitted. \labsw{}{(we have explicitly limited this for this lab by setting \ishell{queue_size=2} for Vesta)} No input files are written or job submissions made until the ``\ishell{run_project}'' function is reached:
\begin{lstlisting}[style=Python]
run_project(sims)
\end{lstlisting}
\noindent
All of the simulations objects have been collected into a list (\ishell{sims}) for submission.
All of the simulation objects have been collected into a list (\ishell{sims}) for submission.
As written, \ishell{O_dimer.py} will only perform calculations at the equilibrium separation distance of 1.2074 Angstrom, since the list of scaling factors (representing stretching or compressing the dimer) only contains one value (\ishell{scales = [1.00]}). Modify the file now to perform DMC calculations across a range of separation distances with each DMC run using the Jastrow factor optimized at the equilibrium separation distance. Specifically, you will want to change the list of scaling factors to include both compression (\ishell{scale<1.0}) and stretch (\ishell{scale>1.0}):
As written, \ishell{O_dimer.py} will perform calculations only at the equilibrium separation distance of 1.2074 {\AA} since the list of scaling factors (representing stretching or compressing the dimer) contains only one value (\ishell{scales = [1.00]}). Modify the file now to perform DMC calculations across a range of separation distances with each DMC run using the Jastrow factor optimized at the equilibrium separation distance. Specifically, you will want to change the list of scaling factors to include both compression (\ishell{scale<1.0}) and stretch (\ishell{scale>1.0}):
\begin{lstlisting}[style=Python]
scales = [1.00,0.90,0.95,1.05,1.10]
\end{lstlisting}
@ -1117,7 +1117,7 @@ Project starting
\noindent
In this case, five simulation ``cascades'' (workflows) have been identified, each one starting and ending with ``\ishell{dft}'' and ``\ishell{qmc}'' runs, respectively. The six status flags (\ishell{setup}, \ishell{sent_files}, \ishell{submitted}, \ishell{finished}, \ishell{got_output}, \ishell{analyzed}) each show \ishell{0}, indicating that no work has been done yet.
Now change ``\ishell{status_only}'' back to \ishell{0}, set ``\ishell{generate_only}'' to \ishell{1}, and run \ishell{O_dimer.py} again. This will perform a dry-run of all simulations. The dry-run should finish in about 20 seconds:
Now change ``\ishell{status_only}'' back to \ishell{0}, set ``\ishell{generate_only}'' to \ishell{1}, and run \ishell{O_dimer.py} again. This will perform a dry run of all simulations. The dry run should finish in about 20 seconds:
\ifws
\begin{shade}
Project starting
@ -1287,7 +1287,7 @@ Project finished
\fi
\noindent
Nexus polls the simulation status every 3 seconds and sleeps in between. The ``scale\_*'' directories should now contain several files:
Nexus polls the simulation status every 3 seconds and sleeps in between. The ``\texttt{scale\_*}'' directories should now contain several files:
\ifws
\begin{shade}
scale_1.0
@ -1367,9 +1367,9 @@ Project starting
Project finished
\end{shade}
\noindent
This way one can continue to add to the \ishell{O_dimer.py} file (\emph{e.g.} adding more separation distances) without worrying about duplicate job submissions.
This way you can continue to add to the \ishell{O_dimer.py} file (e.g., adding more separation distances) without worrying about duplicate job submissions.
Let's actually submit the jobs in the dimer workflow now. Reset the state of the simulations by removing the \ishell{sim.p} files (``\ishell{rm ./scale*/sim*/sim.p}''), set ``\ishell{generate_only}'' to \ishell{0}, and rerun \ishell{O_dimer.py}. It should take about 20 minutes for all the jobs to complete. You may wish to open another terminal to monitor the progress of the individual jobs while the current terminal runs \ishell{O_dimer.py} in the foreground. You can begin the first exercise below once the optimization job completes.
Now submit the jobs in the dimer workflow. Reset the state of the simulations by removing the \ishell{sim.p} files (``\ishell{rm ./scale*/sim*/sim.p}''), set ``\ishell{generate_only}'' to \ishell{0}, and rerun \ishell{O_dimer.py}. It should take about 20 minutes for all the jobs to complete. You may wish to open another terminal to monitor the progress of the individual jobs while the current terminal runs \ishell{O_dimer.py} in the foreground. You can begin the following first exercise once the optimization job completes.
\vspace{3cm}
\begin{flushleft}
@ -1380,26 +1380,26 @@ Let's actually submit the jobs in the dimer workflow now. Reset the state of th
\item{Evaluate the traces of the local energy and the DMC walker population for each separation distance with the \ishell{qmca} tool. Are there any anomalies in the runs? Is the acceptance ratio reasonable? Is the wavefunction of similar quality across all separation distances?}
\item{Use the \ishell{dimer_fit.py} tool located in \ishell{oxygen_dimer} to fit the oxygen dimer binding curve. To get the binding energy of the dimer, we will need the DMC energy of the atom. Before performing the fit, answer: What DMC timestep should be used for the oxygen atom results? The tool accepts three arguments (``\ishell{./dimer_fit.py P N E Eerr}''}), \ishell{P} is the prefix of the DMC input files (should be ``\ishell{qmc}'' at this point), \ishell{N} is the order of the fit (use 2 to start), \ishell{E} and \ishell{Eerr} are your DMC total energy and error bar, respectively for the oxygen atom (in eV). A plot of the dimer data will be displayed and text output will show the DMC equilibrium bond length and binding energy as well as experimental values. How accurately does your fit to the DMC data reproduce the experimental values? What factors affect the accuracy of your results?
\item{Use the \ishell{dimer_fit.py} tool located in \ishell{oxygen_dimer} to fit the oxygen dimer binding curve. To get the binding energy of the dimer, we will need the DMC energy of the atom. Before performing the fit, answer: What DMC time step should be used for the oxygen atom results? The tool accepts three arguments (``\ishell{./dimer_fit.py P N E Eerr}''}), \ishell{P} is the prefix of the DMC input files (should be ``\ishell{qmc}'' at this point), \ishell{N} is the order of the fit (use 2 to start), \ishell{E} and \ishell{Eerr} are your DMC total energy and error bar, respectively, for the oxygen atom (in electron volts). A plot of the dimer data will be displayed, and text output will show the DMC equilibrium bond length and binding energy as well as experimental values. How accurately does your fit to the DMC data reproduce the experimental values? What factors affect the accuracy of your results?
\item{Refit your data with a fourth-order polynomial. How do your predictions change with a fourth-order fit? Is a fourth-order fit appropriate for the available data?}
\item{Add new ``\ishell{scale}'' values to the list in \ishell{O_dimer.py} that interpolate between the original set (e.g. expand to \ishell{[1.00,0.90,0.925,0.95,0.975,1.025,1.05,1.075,1.10]}). Perform the DMC calculations and redo the fits. How accurately does your fit to the DMC data reproduce the experimental values? Should this pseudopotential be used in production calculations?}
\item{Add new ``\ishell{scale}'' values to the list in \ishell{O_dimer.py} that interpolate between the original set (e.g., expand to \ishell{[1.00,0.90,0.925,0.95,0.975,1.025,1.05,1.075,1.10]}). Perform the DMC calculations and redo the fits. How accurately does your fit to the DMC data reproduce the experimental values? Should this pseudopotential be used in production calculations?}
\item{(Optional) Perform optimization runs at the extremal separation distances corresponding to \ishell{scale=[0.90,1.10]}}. Are the individually optimized wavefunctions of significantly better quality than the one imported from \ishell{scale=1.00}? Why? What form of Jastrow factor might give an even better improvement?
\item{(Optional) Perform optimization runs at the extreme separation distances corresponding to \ishell{scale=[0.90,1.10]}}. Are the individually optimized wavefunctions of significantly better quality than the one imported from \ishell{scale=1.00}? Why? What form of Jastrow factor might give an even better improvement?
\end{enumerate}
\section{(Optional) Running your system with QMCPACK}\label{sec:your_system}
This section covers a fairly simple route to get started on QMC calculations of an arbitrary system of interest using the Nexus workflow management system to setup input files and optionally perform the runs. The example provided in this section uses QM Espresso (PWSCF) to generate the orbitals forming the Slater determinant part of the trial wavefunction. PWSCF is a natural choice for solid state systems and it can be used for surface/slab and molecular systems as well, albeit at the price of describing additional vacuum space with plane waves.
This section covers a fairly simple route to get started on QMC calculations of an arbitrary system of interest using the Nexus workflow management system to set up input files and optionally perform the runs. The example provided in this section uses QE (PWSCF) to generate the orbitals forming the Slater determinant part of the trial wavefunction. PWSCF is a natural choice for solid-state systems, and it can be used for surface/slab and molecular systems as well, albeit at the price of describing additional vacuum space with plane waves.
To start out with, you will need pseudopotentials (PP's) for each element in your system in both the UPF (PWSCF) and FSATOM/XML (QMCPACK) formats. A good place to start is the Burkatzki-Filippi-Dolg (BFD) pseudopotential database \newline (\href{http://www.burkatzki.com/pseudos/index.2.html}{http://www.burkatzki.com/pseudos/index.2.html}), which we have already used in our study of the oxygen atom. The database does not contain PP's for the 4th and 5th row transition metals or any of the lanthanides or actinides. If you need a PP that is not in the BFD database, you may need to generate and test one manually (\emph{e.g.} with OPIUM, \href{http://opium.sourceforge.net/}{http://opium.sourceforge.net/}). Otherwise, use \ishell{ppconvert} as outlined in section~\ref{sec:lqb_pseudo} to obtain PP's in the formats used by PWSCF and QMCPACK. Enter the \ishell{your_system} lab directory and place the converted PP's in \ishell{your_system/pseudopotentials}.
To start out, you will need PPs for each element in your system in both the UPF (PWSCF) and FSATOM/XML (QMCPACK) formats. A good place to start is the BFD pseudopotential database \newline (\href{http://www.burkatzki.com/pseudos/index.2.html}{http://www.burkatzki.com/pseudos/index.2.html}), which we have already used in our study of the oxygen atom. The database does not contain PPs for the fourth and fifth row transition metals or any of the lanthanides or actinides. If you need a PP that is not in the BFD database, you may need to generate and test one manually (e.g., with OPIUM, \href{http://opium.sourceforge.net/}{http://opium.sourceforge.net/}). Otherwise, use \ishell{ppconvert} as outlined in Section \ref{sec:lqb_pseudo} to obtain PPs in the formats used by PWSCF and QMCPACK. Enter the \ishell{your_system} lab directory and place the converted PPs in \ishell{your_system/pseudopotentials}.
Before performing production calculations (more than just the initial setup in this section) be sure to converge the plane wave energy cutoff in PWSCF as these PP's can be rather hard, sometimes requiring cutoffs in excess of 300 Ry. Depending on the system under study, the amount of memory required to represent the orbitals (QMCPACK uses 3D B-splines) becomes prohibitive and one may be forced to search for softer PP's.
Before performing production calculations (more than just the initial setup in this section), be sure to converge the plane-wave energy cutoff in PWSCF as these PPs can be rather hard, sometimes requiring cutoffs in excess of 300 Ry. Depending on the system under study, the amount of memory required to represent the orbitals (QMCPACK uses 3D B-splines) can become prohibitive, forcing you to search for softer PPs.
Beyond pseudopotentials, all that is required to get started are the atomic positions and the dimensions/shape of the simulation cell. The Nexus file \ishell{example.py} illustrates how to setup PWSCF and QMCPACK input files by providing minimal information regarding the physical system (an 8-atom cubic cell of diamond in the example). Most of the contents should be familiar from your experience with the automated calculations of the oxygen dimer binding curve in section~\ref{sec:dimer_automation} (if you've skipped ahead you may want to skim that section for relevant information). The most important change is the expanded description of the physical system:
Beyond PPs, all that is required to get started are the atomic positions and the dimensions/shape of the simulation cell. The Nexus file \ishell{example.py} illustrates how to set up PWSCF and QMCPACK input files by providing minimal information regarding the physical system (an 8-atom cubic cell of diamond in the example). Most of the contents should be familiar from your experience with the automated calculations of the oxygen dimer binding curve in Section \ref{sec:dimer_automation} (if you have skipped ahead you may want to skim that section for relevant information). The most important change is the expanded description of the physical system:
\begin{lstlisting}[style=Python]
# details of your physical system (diamond conventional cell below)
@ -1442,7 +1442,7 @@ my_bconds = 'ppp' # ppp/nnn for periodic/open BC's in QMC
# if nnn, center atoms about (a1+a2+a3)/2
\end{lstlisting}
If you have a system you would like to try with QMC, make a copy of \ishell{example.py} and fill in the relevant information about the pseudopotentials, simulation cell axes, and atomic species/positions. Otherwise, you can proceed with \ishell{example.py} as it is.
If you have a system you would like to try with QMC, make a copy of \ishell{example.py} and fill in the relevant information about the PPs, simulation cell axes, and atomic species/positions. Otherwise, you can proceed with \ishell{example.py} as it is.
%The other new aspects are two additional compute jobs to generate the orbitals with PWSCF and convert them into the ESHDF format with \ishell{pw2qmcpack.x}:
@ -1495,7 +1495,7 @@ qsub --mode script --env BG_SHAREDMEMSIZE=32 scf.qsub.in
qsub --mode script --env BG_SHAREDMEMSIZE=32 p2q.qsub.in
\end{shade}
\fi
Once the conversion process has finished the orbitals should be located in the file \newline \ishell{diamond_vmc/pwscf_output/pwscf.pwscf.h5}. Open \ishell{diamond_vmc/vmc.in.xml} and replace ``\ishell{MISSING.h5}'' with ``\ishell{./pwscf_output/pwscf.pwscf.h5}''. Next submit the VMC run:
Once the conversion process has finished, the orbitals should be located in the file \newline \ishell{diamond_vmc/pwscf_output/pwscf.pwscf.h5}. Open \ishell{diamond_vmc/vmc.in.xml} and replace ``\ishell{MISSING.h5}'' with ``\ishell{./pwscf_output/pwscf.pwscf.h5}''. Next submit the VMC run:
\ifws
\begin{shade}
mpirun -np 4 qmcpack vmc.in.xml>&vmc.out&
@ -1505,7 +1505,7 @@ mpirun -np 4 qmcpack vmc.in.xml>&vmc.out&
qsub --mode script --env BG_SHAREDMEMSIZE=32 vmc.qsub.in
\end{shade}
\fi
Note: If your system is large, the above process may not complete within the time frame of this lab. Working with a stripped down (but relevant) example is a good idea for exploratory runs.
Note: If your system is large, the preceding process may not complete within the time frame of this lab. Working with a stripped down (but relevant) example is a good idea for exploratory runs.
Once the runs have finished, you may want to begin exploring Jastrow optimization and DMC for your system. Example calculations are provided at the end of \ishell{example.py} in the commented out text.
@ -1524,7 +1524,7 @@ Once the runs have finished, you may want to begin exploring Jastrow optimizatio
%\hide{
\section{Appendix A: Basic Python constructs\label{app:python_basics}}
Basic Python data types (\ishell{int}, \ishell{float}, \ishell{str}, \ishell{tuple}, \ishell{list}, \ishell{array}, \ishell{dict}, \ishell{obj}) and programming constructs (\ishell{if} statements, \ishell{for} loops, functions w/ keyword arguments) are briefly overviewed below. All examples can be executed interactively in Python. To do this, type ``\ishell{python}'' at the command line and paste any of the shaded text below at the ``\ishell{>>>}'' prompt. For more information about effective use of Python, consult the detailed online documentation: \href{https://docs.python.org/2/}{https://docs.python.org/2/}.
Basic Python data types (\ishell{int}, \ishell{float}, \ishell{str}, \ishell{tuple}, \ishell{list}, \ishell{array}, \ishell{dict}, \ishell{obj}) and programming constructs (\ishell{if} statements, \ishell{for} loops, functions w/ keyword arguments) are briefly overviewed in the following. All examples can be executed interactively in Python. To do this, type ``\ishell{python}'' at the command line and paste any of the shaded text below at the ``\ishell{>>>}'' prompt. For more information about effective use of Python, consult the detailed online documentation: \href{https://docs.python.org/2/}{https://docs.python.org/2/}.
\subsubsection{Intrinsic types: \ishell{int, float, str}}
\begin{shade}
@ -1601,7 +1601,7 @@ o.c = 7 # add/assign element
o.set(c=7,d=8) # add/assign multiple elements
\end{shade}
An important feature of Python to be aware of is that assignment is most often by reference, \emph{i.e.} new values are not always created. This point is illustrated below with an \ishell{obj} instance, but it also holds for \ishell{list}, \ishell{array}, \ishell{dict}, and others.
An important feature of Python to be aware of is that assignment is most often by reference, that is, new values are not always created. This point is illustrated with an \ishell{obj} instance in the following example, but it also holds for \ishell{list}, \ishell{array}, \ishell{dict}, and others.
\begin{shade}
>>> o = obj(a=5,b=6)
>>>

265
manual/lab_qmc_statistics.tex
View File

@ -1,28 +1,27 @@
\chapter{Lab 1: Monte Carlo Statistical Analysis}
\chapter{Lab 1: MC Statistical Analysis}
\label{chap:lab_qmc_statistics}
\section{Topics covered in this Lab}
\section{Topics covered in this lab}
This lab focuses on the basics of analyzing data from Monte Carlo (MC)
This lab focuses on the basics of analyzing data from MC
calculations. In this lab, participants will use data from
VMC calculations of a simple one-electron system with an analytically soluble
system (the ground state of the hydrogen atom) to understand how to interpret a
MC situation. Most of these analyses will also carry over to diffusion Monte
Carlo (DMC) simulations. Topics covered include:
VMC calculations of a simple 1-electron system with an analytically soluble
system (the ground state of the hydrogen atom) to understand how to interpret an
MC situation. Most of these analyses will also carry over to DMC simulations. Topics covered include:
\begin{itemize}
\item{averaging Monte Carlo variables}
\item{the statisical error bar of mean values}
\item{effects of autocorrelation and variance on the error bar}
\item{the relationship between Monte Carlo timestep and autocorrelation}
\item{the use of blocking to reduce autocorrelation}
\item{the significance of the acceptance ratio}
\item{the significance of the sample size}
\item{how to determine whether a Monte Carlo run was successful}
\item{the relationship between wavefunction quality and variance}
\item{gauging the efficiency of Monte Carlo runs}
\item{the cost of scaling up to larger system sizes}
\item{Averaging MC variables}
\item{The statisical error bar of mean values}
\item{The effects of autocorrelation and variance on the error bar}
\item{The relationship between MC time step and autocorrelation}
\item{The use of blocking to reduce autocorrelation}
\item{The significance of the acceptance ratio}
\item{The significance of the sample size}
\item{How to determine whether an MC run was successful}
\item{The relationship between wavefunction quality and variance}
\item{Gauging the efficiency of MC runs}
\item{The cost of scaling up to larger system sizes}
\end{itemize}
@ -141,11 +140,11 @@ labs/lab1_qmc_statistics/
\section{Atomic units}
QMCPACK operates in Hartree atomic units to reduce the
QMCPACK operates in Ha atomic units to reduce the
number of factors in the Schr\"odinger equation. Thus, the unit of length is
the bohr (5.291772 $\times 10^{-11}$ m = 0.529177 \AA); the unit of energy is
the hartree (4.359744 $\times 10^{-18}$ J = 27.211385 eV). The energy of the
ground state of the hydrogen atom in these units is -0.5 hartrees.
the Ha (4.359744 $\times 10^{-18}$ J = 27.211385 eV). The energy of the
ground state of the hydrogen atom in these units is -0.5 Ha.
%\section{Monte Carlo data analysis:\newline average, error bars, variance}
@ -153,18 +152,17 @@ ground state of the hydrogen atom in these units is -0.5 hartrees.
\section{Reviewing statistics}
\label{sec:review}
We will practice taking the average (mean) and standard deviation of some Monte
Carlo data by hand to review the basic definitions.
We will practice taking the average (mean) and standard deviation of some MC data by hand to review the basic definitions.
Enter Python's command line by typing \textbf{python [Enter]}.
You will see a prompt ``\textgreater\textgreater\textgreater''.
Enter Python's command line by typing \texttt{python [Enter]}.
You will see a prompt ``\textgreater\textgreater\textgreater.''
The mean of a data set is given by:
The mean of a dataset is given by:
\begin{align}
\overline{x} = \frac{1}{N}\sum_{i=1}^{N} x_i
\overline{x} = \frac{1}{N}\sum_{i=1}^{N} x_i\:.
\end{align}
To calculate the average of five local energies from a MC calculation of the
To calculate the average of five local energies from an MC calculation of the
ground state of an electron in the hydrogen atom, input (truncate at the
thousandths place if you cannot copy and paste; script versions are also
available in the \texttt{average} directory):
@ -179,7 +177,7 @@ available in the \texttt{average} directory):
)/5.
\end{lstlisting}
Then, press \textbf{[Enter]} to get:
Then, press \texttt{[Enter]} to get:
\begin{shade}
>>> ((-0.45298911858) + (-0.45481953564) + (-0.48066105923) +
@ -191,7 +189,7 @@ To understand the significance of the mean, we also need the standard deviation
around the mean of the data (also called the error bar), given by:
\begin{align}
\sigma = \sqrt{\frac{1}{N(N-1)}\sum_{i=1}^{N} ({x_i} - \overline{x})^2}
\sigma = \sqrt{\frac{1}{N(N-1)}\sum_{i=1}^{N} ({x_i} - \overline{x})^2}\:.
\end{align}
To calculate the standard deviation around the mean (-0.464736835668) of these
@ -207,7 +205,7 @@ five data points, put in:
)**0.5
\end{lstlisting}
Then, press \textbf{[Enter]} to get:
Then, press \texttt{[Enter]} to get:
\begin{shade}
>>> ( (1./(5.*(5.-1.))) * ( (-0.45298911858-(-0.464736835668))**2 +
@ -217,14 +215,14 @@ Then, press \textbf{[Enter]} to get:
0.0053303187464332066
\end{shade}
Thus, we might report this data as having a value -0.465 +/- 0.005 hartrees.
Thus, we might report this data as having a value -0.465 +/- 0.005 Ha.
This calculation of the standard deviation assumes that the average for this
data is fixed, but we may continually add Monte Carlo samples to the data so it
data is fixed, but we can continually add MC samples to the data, so it
is better to use an estimate of the error bar that does not rely on the overall
average. Such an estimate is given by:
\begin{align}
\tilde{\sigma} = \sqrt{\frac{1}{N-1}\sum_{i=1}^{N} \left[{(x^2)}_i - ({x_i})^2\right]}
\tilde{\sigma} = \sqrt{\frac{1}{N-1}\sum_{i=1}^{N} \left[{(x^2)}_i - ({x_i})^2\right]}\:.
\end{align}
To calculate the standard deviation with this formula, input the following,
@ -241,7 +239,7 @@ corresponding local energy:
)**0.5
\end{lstlisting}
and press \textbf{[Enter]} to get:
and press \texttt{[Enter]} to get:
\begin{shade}
>>> ((1./(5.-1.))*((0.60984565298-(-0.45298911858)**2)+
@ -252,20 +250,20 @@ and press \textbf{[Enter]} to get:
\end{shade}
This much larger standard deviation, acknowledging that the mean of this small
data set is not the average in the limit of infinite sampling more accurately,
reports the value of the local energy as -0.5 +/- 0.8 hartrees.
data set is not the average in the limit of infinite sampling, more accurately
reports the value of the local energy as -0.5 +/- 0.8 Ha.
Type \textbf{quit()} and press \textbf{[Enter]} to exit the Python command line.
Type \texttt{quit()} and press \texttt{[Enter]} to exit the Python command line.
\section{Inspecting Monte Carlo data}
\section{Inspecting MC data}
\label{sec:inspect_data}
QMCPACK outputs data from MC calculations into files ending in scalar.dat.
Several quantities are calculated and written for each block of Monte Carlo
QMCPACK outputs data from MC calculations into files ending in \texttt{scalar.dat}.
Several quantities are calculated and written for each block of MC
steps in successive columns to the right of the step index.
Change directories to \texttt{atom}, and open the file ending in
scalar.dat with a text editor (e.g., \textbf{vi *.scalar.dat} or \textbf{emacs
\texttt{scalar.dat} with a text editor (e.g., \textbf{vi *.scalar.dat} or \textbf{emacs
*.scalar.dat}. If possible, adjust the terminal so that lines do not wrap.
The data will begin as follows (broken into three groups to fit on this page):
@ -287,14 +285,13 @@ The first line begins with a \#, indicating that this line does not contain MC
data but rather the labels of the columns. After a blank line, the remaining
lines consist of the MC data. The first column, labeled index, is an integer
indicating which block of MC data is on that line. The second column contains
the quantity usually of greatest interest from the simulation, the local
energy. Since this simulation did not use the exact ground state wave
function, it does not produce -0.5 hartrees as the local energy although the
the quantity usually of greatest interest from the simulation: the local
energy. Since this simulation did not use the exact ground state wavefunction, it does not produce -0.5 Ha as the local energy although the
value lies within about 10\%. The value of the local energy fluctuates from
block to block and the closer the trial wave function is to the ground state,
block to block, and the closer the trial wavefunction is to the ground state
the smaller these fluctuations will be. The next column contains an important
ingredient in estimating the error in the MC average--the square of the local
energy--found by evaluating the square of the Hamiltonian.
ingredient in estimating the error in the MC average---the square of the local
energy---found by evaluating the square of the Hamiltonian.
\begin{shade}
... Kinetic Coulomb BlockWeight ...
@ -317,7 +314,7 @@ only the potential energy coming from its interaction with the nucleus. In
many-electron simulations, the local potential energy contains contributions
from the electron-electron Coulomb interactions and the nuclear potential or
pseudopotential. The fifth column contains the local kinetic energy value for
each MC block, obtained from the Laplacian of the wave function. The sixth
each MC block, obtained from the Laplacian of the wavefunction. The sixth
column shows the local Coulomb interaction energy. The seventh column displays
the weight each line of data has in the average (the weights are identical in
this simulation).
@ -338,7 +335,7 @@ this simulation).
The eighth column shows the CPU time (in seconds) to calculate the data in that
line. The ninth column from the left contains the acceptance ratio (1 being
full acceptance) for Monte Carlo steps in that line's data. Other than the
full acceptance) for MC steps in that line's data. Other than the
block weight, all quantities vary from line to line.
Exit the text editor (\textbf{[Esc] :q! [Enter]} in vi, \textbf{[Ctrl]-x [Ctrl]-c} in
@ -347,11 +344,11 @@ emacs).
\section{Averaging quantities in the MC data}
\label{sec:averaging}
QMCPACK includes the qmca Python tool to average quantities in the scalar.dat file (and
also the dmc.dat file of DMC simulations). Without any flags, qmca will output
the average of each column with a quantity in the scalar.dat file as follows.
QMCPACK includes the qmca Python tool to average quantities in the \texttt{scalar.dat} file (and
also the \texttt{dmc.dat} file of DMC simulations). Without any flags, qmca will output
the average of each column with a quantity in the \texttt{scalar.dat} file as follows.
Execute qmca by \textbf{qmca *.scalar.dat}, which for this data outputs:
Execute qmca by \texttt{qmca *.scalar.dat}, which for this data outputs:
\begin{shade}
@ -369,9 +366,9 @@ Efficiency = 0.00000000 +/- 0.00000000
\end{shade}
After one blank, qmca prints the title of the subsequent data, gleaned from the
data file name. In this case, H.s000.scalar.dat became ``H series 0''.
Everything before the first ``.s'' will be interpreted as the title, and the
number between ``.s'' and the next ``.'' will be interpreted as the series
data file name. In this case, \texttt{H.s000.scalar.dat} became ``H series 0.''
Everything before the first ``\texttt{.s}'' will be interpreted as the title, and the
number between ``\texttt{.s}'' and the next ``.'' will be interpreted as the series
number.
The first column under the title is the name of each quantity qmca averaged.
@ -380,18 +377,18 @@ quantity of that line, and the column to the right of the plus-slash-minus is
the statistical error bar on the quantity. All quantities calculated from MC
simulations have and must be reported with a statistical error bar!
Two new quantities not present in the scalar.dat file are computed by qmca from
the data--variance and efficiency. We will look at these later in this lab.
Two new quantities not present in the \texttt{scalar.dat} file are computed by qmca from
the data---variance and efficiency. We will look at these later in this lab.
To view only one value, \textbf{qmca} takes the \textbf{-q (quantity)} flag.
For example, the output of \textbf{qmca -q LocalEnergy *.scalar.dat} in this
For example, the output of \texttt{qmca -q LocalEnergy *.scalar.dat} in this
directory produces a single line of output:
\begin{shade}
H series 0 LocalEnergy = -0.454460 +/- 0.000568
\end{shade}
Type \textbf{qmca --help} to see the list of all quantities and their
Type \texttt{qmca --help} to see the list of all quantities and their
abbreviations.
\section{Evaluating MC simulation quality}
@ -401,15 +398,15 @@ it went. Besides the deviation of the average from an expected value (if there
is one), the stability of the simulation in its sampling, the autocorrelation
between MC steps, the value of the acceptance ratio (accepted steps over total
proposed steps), and the variance in the local energy all indicate the quality
of a MC simulation. We will look at these one by one.
of an MC simulation. We will look at these one by one.
\subsection{Tracing MC quantities}
Visualizing the evolution of MC quantities over the course of the simulation by
a \textit{trace} offers a quick picture of whether the random walk had expected
a \textit{trace} offers a quick picture of whether the random walk had the expected
behavior. qmca plots traces with the -t flag.
Type \textbf{qmca -q e -t H.s000.scalar.dat}, which produces a graph of the
Type \texttt{qmca -q e -t H.s000.scalar.dat}, which produces a graph of the
trace of the local energy:
\FloatBarrier
@ -427,13 +424,13 @@ The solid black line connects the values of the local energy at each MC block
red line. One standard deviation above and below the average are marked with
horizontal, dashed red lines.
The trace of this run is largely centered around the average with no
large-scale oscillations or major shifts, indicating a good quality MC run.
The trace of this run is largely centered on the average with no
large-scale oscillations or major shifts, indicating a good-quality MC run.
Try tracing the kinetic and potential energies, seeing that their behavior is
comparable to the total local energy.
comparable with the total local energy.
Change to directory \texttt{problematic} and type \textbf{qmca -q e -t
Change to directory \texttt{problematic} and type \texttt{qmca -q e -t
H.s000.scalar.dat} to produce this graph:
\FloatBarrier
@ -446,21 +443,21 @@ H.s000.scalar.dat} to produce this graph:
%\includegraphics[scale=0.5]{E_L_H_B-splines.png}
Here, the local energy samples cluster around the expected -0.5 hartrees for the
Here, the local energy samples cluster around the expected -0.5 Ha for the
first 150 samples or so and then begin to oscillate more wildly and increase
erratically toward 0, indicating a poor quality MC run.
erratically toward 0, indicating a poor-quality MC run.
Again, trace the kinetic and potential energies in this run and see how their
behavior compares to the total local energy.
behavior compares with the total local energy.
\subsection{Blocking away autocorrelation}
\textit{Autocorrelation} occurs when a given MC step biases subsequent MC
steps, leading to samples that are not statistically independent. We must take
this autocorrelation into account in order to obtain accurate statistics. qmca
this autocorrelation into account to obtain accurate statistics. qmca
outputs autocorrelation when given the {-}{-}sac flag.
Change to directory \texttt{autocorrelation} and type \textbf{qmca -q e
Change to directory \texttt{autocorrelation} and type \texttt{qmca -q e
{-}{-}sac H.s000.scalar.dat}.
\begin{shade}
@ -472,8 +469,8 @@ this case).
Proposing too small a step in configuration space, the MC \textit{time step},
can lead to autocorrelation since the new samples will be in the neighborhood
of previous samples. Type \textbf{grep timestep H.xml} to see the varying time
step values in this QMCPACK input file (H.xml):
of previous samples. Type \texttt{grep timestep H.xml} to see the varying time
step values in this QMCPACK input file (\texttt{H.xml}):
\begin{shade}
<parameter name="timestep">10</parameter>
@ -496,8 +493,8 @@ step values in this QMCPACK input file (H.xml):
Generally, as the time step decreases, the autocorrelation will increase
(caveat: very large time steps will also have increasing autocorrelation). To
see this, type \textbf{qmca -q e {-}{-}sac *.scalar.dat} to see the energies
and autocorrelation times, then plot with gnuplot by inputting \textbf{gnuplot
see this, type \texttt{qmca -q e {-}{-}sac *.scalar.dat} to see the energies
and autocorrelation times, then plot with gnuplot by inputting \texttt{gnuplot
H.plt}:
\FloatBarrier
@ -512,25 +509,25 @@ H.plt}:
The error bar also increases with the autocorrelation.
Press \textbf{q [Enter]} to quit gnuplot.
Press \texttt{q [Enter]} to quit gnuplot.
To get around the bias of autocorrelation, we group the MC steps into blocks,
take the average of the data in the steps of each block, and then finally
average the averages in all the blocks. QMCPACK outputs the block averages as
each line in the scalar.dat file. (For DMC simulations, in addition to the
scalar.dat, QMCPACK outputs the quantities at each step to the dmc.dat file,
each line in the \texttt{scalar.dat} file. (For DMC simulations, in addition to the
\texttt{scalar.dat}, QMCPACK outputs the quantities at each step to the \texttt{dmc.dat} file,
which permits reblocking the data differently from the specification in the
input file.)
Change directories to \texttt{blocking}. Here we look at the time step of the
last data set in the \texttt{autocorrelation} directory. Verify this by typing
\textbf{grep timestep H.xml} to see that all values are set to 0.001. Now to
see how we will vary the blocking, type \textbf{grep -A1 blocks H.xml}. The
last dataset in the \texttt{autocorrelation} directory. Verify this by typing
\texttt{grep timestep H.xml} to see that all values are set to 0.001. Now to
see how we will vary the blocking, type \texttt{grep -A1 blocks H.xml}. The
parameter ``steps'' indicates the number of steps per block, and the parameter
``blocks'' gives the number of blocks. For this comparison, the total number
of MC steps (equal to the product of ``steps'' and ``blocks'') is fixed at
50000. Now check the effect of blocking on autocorrelation--type \textbf{qmca
-q e {-}{-}sac *scalar.dat} to see the data and \textbf{gnuplot H.plt} to
50,000. Now check the effect of blocking on autocorrelation---type \texttt{qmca
-q e {-}{-}sac *scalar.dat} to see the data and \texttt{gnuplot H.plt} to
visualize the data:
%\begin{shaded}
@ -569,7 +566,7 @@ time. The larger number of blocks over which to average at small
step-per-block number masks the corresponding increase in error bar with
increasing autocorrelation.
Press \textbf{q [Enter]} to quit gnuplot.
Press \texttt{q [Enter]} to quit gnuplot.
\subsection{Balancing autocorrelation and acceptance ratio}
@ -579,8 +576,8 @@ probability distribution is similar and thus more likely to result in an
accepted move. Keeping the acceptance ratio high means the algorithm is
efficiently exploring configuration space and not sticking at particular
configurations. Return to the \ishell{autocorrelation} directory. Refresh your
memory on the time steps in this set of simulations by \textbf{grep timestep
H.xml}. Then, type \textbf{qmca -q ar *scalar.dat} to see the acceptance ratio
memory on the time steps in this set of simulations by \texttt{grep timestep
H.xml}. Then, type \texttt{qmca -q ar *scalar.dat} to see the acceptance ratio
as it varies with decreasing time step:
\begin{shade}
@ -607,29 +604,28 @@ By series 8 (time step = 0.02), the acceptance ratio is in excess of 99\%.
Considering the increase in autocorrelation and subsequent increase in error
bar as time step decreases, it is important to choose a time step that trades
off appropriately between acceptance ratio and autocorrelation. In this
example, a time step of 0.02 occupies a spot where acceptance ratio is high
(99.6\%), and autocorrelation is not appreciably larger than the minimum value
example, a time step of 0.02 occupies a spot where the acceptance ratio is high
(99.6\%) and autocorrelation is not appreciably larger than the minimum value
(1.4 vs. 1.0).
\subsection{Considering variance}
Besides autocorrelation, the dominant contributor to the error bar is the
\textit{variance} in the local energy. The variance measures the fluctuations
around the average local energy, and, as the fluctuations go to zero, the wave
function reaches an exact eigenstate of the Hamiltonian. qmca calculates this
from the local energy and local energy squared columns of the scalar.dat.
around the average local energy, and, as the fluctuations go to zero, the wavefunction reaches an exact eigenstate of the Hamiltonian. qmca calculates this
from the local energy and local energy squared columns of the \texttt{scalar.dat}.
Type \textbf{qmca -q v H.s009.scalar.dat} to calculate the variance on the run
Type \texttt{qmca -q v H.s009.scalar.dat} to calculate the variance on the run
with time step balancing autocorrelation and acceptance ratio:
\begin{shade}
H series 9 Variance = 0.513570 +/- 0.010589
\end{shade}
Just as the total energy doesn't tell us much by itself, neither does the
variance. However, comparing the ratio of the variance to the energy indicates
how the magnitude of the fluctuations compares to the energy itself. Type
\textbf{qmca -q ev H.s009.scalar.dat} to calculate the energy and variance on
Just as the total energy does not tell us much by itself, neither does the
variance. However, comparing the ratio of the variance with the energy indicates
how the magnitude of the fluctuations compares with the energy itself. Type
\texttt{qmca -q ev H.s009.scalar.dat} to calculate the energy and variance on
the run side by side with the ratio:
\begin{shade}
@ -637,23 +633,23 @@ the run side by side with the ratio:
H series 0 -0.454460 +/- 0.000568 0.529496 +/- 0.018445 1.1651
\end{shade}
1.1651 is a very high ratio indicating the square of the fluctuations is on
The very high ration of 1.1651 indicates the square of the fluctuations is on
average larger than the value itself. In the next section, we will approach
ways to improve the variance that subsequent labs will build upon.
ways to improve the variance that subsequent labs will build on.
\section{Reducing statistical error bars}
\subsection{Increasing MC sampling}
Increasing the number of MC samples in a data set reduces the error bar as the
Increasing the number of MC samples in a dataset reduces the error bar as the
inverse of the square root of the number of samples. There are two ways to
increase the number of MC samples in a simulation: running more samples in
parallel and increasing the number of blocks (with fixed number of steps per
increase the number of MC samples in a simulation: (1) running more samples in
parallel and (2) increasing the number of blocks (with fixed number of steps per
block, this increases the total number of MC steps).
To see the effect of the running more samples in parallel, change to the
To see the effect of running more samples in parallel, change to the
directory \ishell{nodes}. The series here increases the number of nodes by
factors of four from 32 to 128 to 512. Type \textbf{qmca -q ev *scalar.dat}
factors of four from 32 to 128 to 512. Type \texttt{qmca -q ev *scalar.dat}
and note the change in the error bar on the local energy as the number of
nodes. Visualize this with \textbf{gnuplot H.plt}:
@ -670,13 +666,13 @@ nodes. Visualize this with \textbf{gnuplot H.plt}:
Increasing the number of blocks, unlike running in parallel, increases the
total CPU time of the simulation.
Press \textbf{q [Enter]} to quit gnuplot.
Press \texttt{q [Enter]} to quit gnuplot.
To see the effect of increasing the block number, change to the directory
\ishell{blocks}. To see how we will vary the number of blocks, type
\textbf{grep -A1 blocks H.xml}. The number of steps remains fixed, thus
\texttt{grep -A1 blocks H.xml}. The number of steps remains fixed, thus
increasing the total number of samples. Visualize the tradeoff by inputting
\textbf{gnuplot H.plt}:
\texttt{gnuplot H.plt}:
\FloatBarrier
\begin{figure}[ht!]
@ -688,22 +684,21 @@ increasing the total number of samples. Visualize the tradeoff by inputting
%\includegraphics[scale=1.0]{nblock_vs_tcpu_energy_H_STO-2G.png}
Press \textbf{q [Enter]} to quit gnuplot.
Press \texttt{q [Enter]} to quit gnuplot.
\subsection{Improving the basis set}
In all of the above examples, we are using the sum of two gaussian functions
(STO-2G) to approximate what should be a simple decaying exponential (STO =
Slater-type orbital) for the wave function of the ground state of the hydrogen
In all of the previous examples, we are using the sum of two Gaussian functions
(STO-2G) to approximate what should be a simple decaying exponential for the wavefunction of the ground state of the hydrogen
atom. The sum of multiple copies of a function varying each copy's width and
amplitude with coefficients is called a \textit{basis set}. As we add gaussians
to the basis set, the approximation improves, the variance goes toward zero and
the energy goes to -0.5 hartrees. In nearly every other case, the exact
amplitude with coefficients is called a \textit{basis set}. As we add Gaussians
to the basis set, the approximation improves, the variance goes toward zero, and
the energy goes to -0.5 Ha. In nearly every other case, the exact
function is unknown, and we add basis functions until the total energy does not
change within some threshold.
Change to the directory \ishell{basis} and look at the total energy and
variance as we change the wave function by typing \textbf{qmca -q ev H\_*}:
variance as we change the wavefunction by typing \texttt{qmca -q ev H\_*}:
\begin{shade}
LocalEnergy Variance ratio
@ -714,12 +709,12 @@ H__exact series 0 -0.500000 +/- 0.000000 0.000000 +/- 0.000000 -0.0000
\end{shade}
qmca also puts out the ratio of the variance to the local energy in a column to
the right of the variance error bar. A typical high quality value for this
ratio is lower than 0.1 or so--none of these few-gaussian wave functions
the right of the variance error bar. A typical high-quality value for this
ratio is lower than 0.1 or so---none of these few-Gaussian wavefunctions
satisfy that rule of thumb.
Use qmca to plot the trace of the local energy, kinetic energy, and potential
energy of H\_\_exact--the total energy is constantly -0.5 hartree even though
energy of H\_\_exact. The total energy is constantly -0.5 Ha even though
the kinetic and potential energies fluctuate from configuration to
configuration.
@ -727,12 +722,12 @@ configuration.
Another route to reducing the variance is the introduction of a Jastrow factor to
account for electron-electron correlation (not the statistical autocorrelation
of Monte Carlo steps but the physical avoidance that electrons have of one another).
of MC steps but the physical avoidance that electrons have of one another).
To do this, we will switch to the hydrogen dimer with the exact ground state
wave function of the atom (STO basis)--this will not be exact for the dimer.
The ground state energy of the hydrogen dimer is -1.174 hartrees.
wavefunction of the atom (STO basis)---this will not be exact for the dimer.
The ground state energy of the hydrogen dimer is -1.174 Ha.
Change directories to \ishell{dimer} and put in \textbf{qmca -q ev *scalar.dat}
Change directories to \ishell{dimer} and put in \texttt{qmca -q ev *scalar.dat}
to see the result of adding a simple, one-parameter Jastrow to the STO basis
for the hydrogen dimer at experimental bond length:
@ -742,9 +737,9 @@ H2_STO___no_jastrow series 0 -0.876548 +/- 0.005313 0.473526 +/- 0.014910
H2_STO_with_jastrow series 0 -0.912763 +/- 0.004470 0.279651 +/- 0.016405
\end{shade}
The energy reduces by 0.044 +/- 0.006 hartrees and the variance by 0.19 +/- 0.02.
The energy reduces by 0.044 +/- 0.006 HA and the variance by 0.19 +/- 0.02.
This is still 20\% above the ground state energy, and subsequent labs will cover how
to improve on this with improved forms of the wave function that capture more
to improve on this with improved forms of the wavefunction that capture more
of the physics.
\section{Scaling to larger numbers of electrons}
@ -753,7 +748,7 @@ of the physics.
The inverse of the product of CPU time and the variance measures the
\textit{efficiency} of an MC calculation. Use qmca to calculate efficiency by
typing \textbf{qmca -q eff *scalar.dat} to see the efficiency of these two
typing \texttt{qmca -q eff *scalar.dat} to see the efficiency of these two
H$_2$ calculations:
\begin{shade}
@ -768,12 +763,12 @@ CPU time to verify this claim).
\subsection{Scaling up}
To see how MC scales with increasing particle number, change directories to
\ishell{size}. Here are the data from runs of increasing number of electrons
\ishell{size}. Here are the data from runs of increasing numbers of electrons
for H, H$_2$, C, CH$_4$, C$_2$, C$_2$H$_4$, (CH$_4$)$_2$, and (C$_2$H$_4$)$_2$
using the STO-6G basis set for the orbitals of the Slater determinant. The file names begin with the number of electrons simulated for those data.
Use \textbf{qmca -q bc *scalar.dat} to see that the CPU time per block
increases with number of electrons in the simulation, then plot the total CPU
Use \texttt{qmca -q bc *scalar.dat} to see that the CPU time per block
increases with the number of electrons in the simulation; then plot the total CPU
time of the simulation by \textbf{gnuplot Nelectron\_tCPU.plt}:
\FloatBarrier
@ -787,13 +782,13 @@ time of the simulation by \textbf{gnuplot Nelectron\_tCPU.plt}:
%\includegraphics[scale=1.0]{nelectron_vs_tcpu_H_C_CH_STO-6G.png}
The green pluses represent the CPU time per block at each electron number.
The red line is a quadratic fit to those data. For a fixed basis set size, we expect the time to scale quadratically up to 1000s of electrons, at which point a cubic scaling term may become dominant. Knowing the scaling allows you to roughly project the calculation time for a larger number of electrons.
The red line is a quadratic fit to those data. For a fixed basis set size, we expect the time to scale quadratically up to 1,000s of electrons, at which point a cubic scaling term may become dominant. Knowing the scaling allows you to roughly project the calculation time for a larger number of electrons.
Press \textbf{q [Enter]} to quit gnuplot.
Press \texttt{q [Enter]} to quit gnuplot.
This isn't the whole story, however. The variance of the energy also increases
This is not the whole story, however. The variance of the energy also increases
with a fixed basis set as the number of particles increases at a faster rate
than the energy decreases. To see this, type \textbf{qmca -q ev *scalar.dat}:
than the energy decreases. To see this, type \texttt{qmca -q ev *scalar.dat}:
\begin{shade}
LocalEnergy Variance
@ -808,5 +803,5 @@ than the energy decreases. To see this, type \textbf{qmca -q ev *scalar.dat}:
\end{shade}
The increase in variance is not uniform, but the general trend is upward with a
fixed wave function form and basis set. Subsequent labs will address how to
improve the wave function in order to keep the variance manageable.
fixed wavefunction form and basis set. Subsequent labs will address how to
improve the wavefunction to keep the variance manageable.

20
manual/pyscf_for_afqmc.tex
View File

@ -5,9 +5,9 @@ The main purpose of pyscf in the group is in the generation of 2-electron integr
For a detailed description of pyscf, please see the documentation by the authors.
Examples of how to generate integrals are:
Examples of how to generate integrals:
\begin{itemize}
\item Run a standard pyscf calculation, e.g. a HF or DFT calculation. Make sure you preserve the chkfile and make sure you store the core hamiltonian on the chkfile. An example of how to do this for a single k-point calculation is found below.
\item Run a standard pyscf calculation, e.g., an HF or DFT calculation. Make sure you preserve the chkfile and make sure you store the core hamiltonian on the chkfile. An example of how to do this for a single k-point calculation follows.
\begin{lstlisting}[style=Python,caption=The following is an example PySCF input file for single k-point calculations.]
import numpy
import h5py
@ -46,7 +46,7 @@ with h5py.File(mf.chkfile) as fh5:
\end{lstlisting}
\item {For a calculation with k-points:
Run a standard pyscf calculation, e.g. a HF or DFT calculation. Make sure you preserve the chkfile and make sure you store the core hamiltonian on the chkfile. An example of how to do this for a single k-point calculation is found below.}
Run a standard pyscf calculation, e.g., an HF or DFT calculation. Make sure you preserve the chkfile and make sure you store the core hamiltonian on the chkfile. An example of how to do this for a single k-point calculation follows.}
\begin{lstlisting}[style=Python,caption=The following is an example PySCF input file for calculations with k-points.]
import numpy
@ -89,7 +89,7 @@ with h5py.File(mf.chkfile) as fh5:
\end{lstlisting}
\end{itemize}
Once the checkpoint file has been created, it is now possible to generate the integral file. The recommended approach is:
Once the checkpoint file has been created, it is possible to generate the integral file. The recommended approach follows:
\begin{lstlisting}[style=Python,caption=The following is an example input file for calculating the integrals.]
from mpi4py import MPI
@ -107,26 +107,26 @@ if rank==0:
integrals_from_chkfile.combine_eri_h5("fcidump", nproc)
\end{lstlisting}
It is also possible to generated the Cholesky decomposed integrals in pyscf directly. This is typically faster and more appropriate. (the phfmol code expects integrals directly, but these can be generated from qmcpack with the Cholesky decomposed ones.) To generate them for a run with kpoints for example, the recommended approach is: (NOTE: for some reason, mpi4py must be imported first!!!)
It is also possible to generated the Cholesky decomposed integrals in pyscf directly. This is typically faster and more appropriate. (The phfmol code expects integrals directly, but these can be generated from QMCPACK with the Cholesky decomposed ones.) To generate them for a run with k-points, for example, the recommended approach follows: (NOTE: for some reason, mpi4py must be imported first!!!)
For calculations with kpoints (those generated with K...), use integrals\_from\_chkfile.eri\_to\_h5\_kpts(...).
Additional arguments to eri\_to\_h5 are:
\begin{itemize}
\item \textbf{cholesky}. Determines whether 2-electron integrals or their cholesky factorization is calculated.
\item \textbf{Cholesky}. Determines whether 2-electron integrals or their Cholesky factorization is calculated.
Default: False
\item \textbf{orthoAO}. If True, generates the integrals in the orthogonalized AO basis. If False, generates the integrals in the MO basis found on the checkpoint file. For UHF calculations, only orthoAO=True is allowed. If set to False, the fock matrix must be stored in the scf dump file.
Default: False
\item \textbf{LINDEP\_CUTOFF}. Cutoff used to define linearly dependent basis functions.
Default: 1e-9
\item \textbf{gtol}. Cutoff applied during writing for 2-electron integrals. If cholesky=True, then this is the cutoff used in the iterative cholesky factorization. In this case, the resulting factorized hamiltonian will have a residual error smaller than the requested cutoff.
\item \textbf{gtol}. Cutoff applied during writing for 2-electron integrals. If Cholesky=True, then this is the cutoff used in the iterative Cholesky factorization. In this case, the resulting factorized Hamiltonian will have a residual error smaller than the requested cutoff.
Default: 1e-6
\item \textbf{wfnName}. Name of the file with the wavefuntion. This is only generated when orthoAO=True.
\item \textbf{wfnName}. Name of the file with the wavefunction. This is generated only when orthoAO=True.
Default: ``wfn.dat''
\item \textbf{wfnPHF}. Name of file with initial guess for phfmol code.
Default: None (no file is generated)
\item \textbf{MaxIntgs}. Maximum number of integrals (or terms in cholesky matrix) per block in the hdf5 file. This controls the size of the hdf5 data sets.
\item \textbf{MaxIntgs}. Maximum number of integrals (or terms in Cholesky matrix) per block in the hdf5 file. This controls the size of the hdf5 datasets.
Default: 2000000
\item \textbf{maxvecs}. Represents the maximum number of Cholesky vectors allowed. The actual maximum number of cholesky vectors in the calculation is set to maxvecsnmo. So a value of 10 leads to an actual cutoff in the number of vectors of 10*nmo, where nmo is the total number of molecular orbitals in the calculation. The calculation will stop at the requested number of vectors even if the tolerance is not reached.
\item \textbf{maxvecs}. Represents the maximum number of Cholesky vectors allowed. The actual maximum number of Cholesky vectors in the calculation is set to maxvecsnmo. So a value of 10 leads to an actual cutoff in the number of vectors of 10*nmo, where nmo is the total number of MOs in the calculation. The calculation will stop at the requested number of vectors even if the tolerance is not reached.
Default: 20
\end{itemize}