.. _lab-qmc-basics: Lab 2: QMC Basics ================= 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: - Converting pseudopotentials into QMCPACK’s FSATOM format - Generating orbitals with QE - Converting orbitals into QMCPACK’s ESHDF format with pw2qmcpack - Optimizing Jastrow factors with QMCPACK - Removing DMC time step errors via extrapolation - Automating QMC workflows with Nexus - Testing pseudopotentials for accuracy Lab outline ----------- #. Download and conversion of oxygen atom pseudopotential #. DMC time step study of the neutral oxygen atom #. DFT orbital generation with QE #. Orbital conversion with #. Optimization of Jastrow correlation factor with QMCPACK #. DMC run with multiple time steps #. DMC time step study of the first ionization potential of oxygen #. Repetition of a-d above for ionized oxygen atom #. Automated DMC calculations of the oxygen dimer binding curve Lab directories and files ------------------------- :: % labs/lab2_qmc_basics/ │ ├── oxygen_atom - oxygen atom calculations │ ├── O.q0.dft.in - Quantum ESPRESSO input for DFT run │ ├── O.q0.p2q.in - pw2qmcpack.x input for orbital conversion run │ ├── O.q0.opt.in.xml - QMCPACK input for Jastrow optimization run │ ├── O.q0.dmc.in.xml - QMCPACK input file for neutral O DMC │ ├── ip_conv.py - tool to fit oxygen IP vs timestep │ └── reference - directory w/ completed runs │ ├── oxygen_dimer - oxygen dimer calculations │ ├── dimer_fit.py - tool to fit dimer binding curve │ ├── O_dimer.py - automation script for dimer calculations │ ├── pseudopotentials - directory for pseudopotentials │ └── reference - directory w/ completed runs │ └── your_system - performing calculations for an arbitrary system (yours) ├── example.py - example nexus file for periodic diamond ├── pseudopotentials - directory containing C pseudopotentials └── reference - directory w/ completed runs .. _lqb-pseudo: Obtaining and converting a pseudopotential for oxygen ----------------------------------------------------- First enter the ``oxygen_atom`` directory: :: cd labs/lab2_qmc_basics/oxygen_atom/ Throughout the rest of the lab, locations are specified with respect to ``labs/lab2_qmc_basics`` (e.g., ``oxygen_atom``). We use a potential from the Burkatzki-Filippi-Dolg pseudopotential database. Although the full database is available in QMCPACK distribution (``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 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 displayed. The desired portion is at the bottom (the last 7 lines). Copy this text into the editor of your choice (e.g., ``emacs`` or ``vi``) and save it as ``O.BFD.gamess`` (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 ``ppconvert`` tool: :: ppconvert --gamess_pot O.BFD.gamess --s_ref "1s(2)2p(4)" \ --p_ref "1s(2)2p(4)" --d_ref "1s(2)2p(4)" --xml O.BFD.xml Observe the notation used to describe the reference valence configuration for this helium-core PP: ``1s(2)2p(4)``. The ``ppconvert`` tool uses the following convention for the valence states: the first $s$ state is labeled ``1s`` (``1s``, ``2s``, ``3s``, :math:`\ldots`), the first :math:`p` state is labeled ``2p`` (``2p``, ``3p``, :math:`\ldots`), and the first :math:`d` state is labeled ``3d`` (``3d``, ``4d``, :math:`\ldots`). Copy the resulting xml file into the ``oxygen_atom`` directory. Note: The command to convert the PP into QE's UPF format is similar (both formats are required): :: ppconvert --gamess_pot O.BFD.gamess --s_ref "1s(2)2p(4)" \ --p_ref "1s(2)2p(4)" --d_ref "1s(2)2p(4)" --log_grid --upf O.BFD.upf For reference, the text of ``O.BFD.gamess`` should be: :: O-QMC GEN 2 1 3 6.00000000 1 9.29793903 55.78763416 3 8.86492204 -38.81978498 2 8.62925665 1 38.41914135 2 8.71924452 The full QMCPACK pseudopotential is also included in ``oxygen_atom/reference/O.BFD.*``. .. _lqb-dft: DFT with QE to obtain the orbital part of the wavefunction ---------------------------------------------------------- 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 ``oxygen_atom/reference`` directory. QE input for the DFT-LDA ground state of the neutral oxygen atom can be found in ``O.q0.dft.in`` and also in :ref:`Listing 58 `. Setting ``wf_collect=.true.`` instructs QE to write the orbitals to disk at the end of the run. Option ``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 :ref:`pw2qmcpack`). Note that the plane-wave energy cutoff has been set to a reasonable value of 300 Ry here (``ecutwfc=300``). This value depends on the pseudopotentials used, and, in general, should be selected by running DFT :math:`\rightarrow` (orbital conversion) :math:`\rightarrow` VMC with increasing energy cutoffs until the lowest VMC total energy and variance is reached. .. code-block:: :caption: QE input file for the neutral oxygen atom (``O.q0.dft.in``) :name: Listing 58 &CONTROL calculation = 'scf' restart_mode = 'from_scratch' prefix = 'O.q0' outdir = './' pseudo_dir = './' disk_io = 'low' wf_collect = .true. / &SYSTEM celldm(1) = 1.0 ibrav = 0 nat = 1 ntyp = 1 nspin = 2 tot_charge = 0 tot_magnetization = 2 input_dft = 'lda' ecutwfc = 300 ecutrho = 1200 nosym = .true. occupations = 'smearing' smearing = 'fermi-dirac' degauss = 0.0001 / &ELECTRONS diagonalization = 'david' mixing_mode = 'plain' mixing_beta = 0.7 conv_thr = 1e-08 electron_maxstep = 1000 / ATOMIC_SPECIES O 15.999 O.BFD.upf ATOMIC_POSITIONS alat O 9.44863067 9.44863161 9.44863255 K_POINTS automatic 1 1 1 0 0 0 CELL_PARAMETERS cubic 18.89726133 0.00000000 0.00000000 0.00000000 18.89726133 0.00000000 0.00000000 0.00000000 18.89726133 Run QE by typing :: mpirun -np 4 pw.x -input O.q0.dft.in >&O.q0.dft.out& The DFT run should take a few minutes to complete. If desired, you can track the progress of the DFT run by typing "``tail -f O.q0.dft.out``." Once finished, you should check the LDA total energy in ``O.q0.dft.out`` by typing "``grep '! ' O.q0.dft.out``." The result should be close to :: ! total energy = -31.57553905 Ry The orbitals have been written in a format native to QE in the ``O.q0.save`` directory. We will convert them into the ESHDF format expected by QMCPACK by using the ``pw2qmcpack.x`` tool. The input for ``pw2qmcpack.x`` can be found in the file ``O.q0.p2q.in`` and also in :ref:`Listing 59 `. .. code-block:: :caption: ``pw2qmcpack.x`` input file for orbital conversion (``O.q0.p2q.in``) :name: Listing 59 &inputpp prefix = 'O.q0' outdir = './' write_psir = .false. / Perform the orbital conversion now by typing the following: :: mpirun -np 1 pw2qmcpack.x&O.q0.p2q.out& Upon completion of the run, a new file should be present containing the orbitals for QMCPACK: ``O.q0.pwscf.h5``. Template XML files for particle (``O.q0.ptcl.xml``) and wavefunction (``O.q0.wfs.xml``) inputs to QMCPACK should also be present. .. _optimization-walkthrough: Optimization with QMCPACK to obtain the correlated part of the wavefunction --------------------------------------------------------------------------- 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: .. math:: :label: eq66 \Psi_{Slater-Jastrow}=e^{-J}\Psi_{Slater} .. math:: :label: eq67 J = \sum_{\sigma\sigma'}\sum_{i p p p 15 with .. code-block:: xml n n n Next we will select Jastrow factors appropriate for an atom. In open boundary conditions, the B-spline Jastrow correlation functions should cut off to zero at some distance away from the atom. Open ``O.q0.wfs.xml`` and add the following cutoffs (``rcut`` in Bohr radii) to the correlation factors: .. code-block:: xml ... ... ... ... These terms correspond to :math:`u_2^{\uparrow\uparrow}/u_2^{\downarrow\downarrow}`, :math:`u_2^{\uparrow\downarrow}`, and :math:`u_1^{\uparrow O}/u_1^{\downarrow O}`, respectively. In each case, the correlation function (:math:`u_*`) is represented by piecewise continuous cubic B-splines. Each correlation function has eight parameters, which are just the values of :math:`u` on a uniformly spaced grid up to ``rcut``. Initially the parameters (``coefficients``) are set to zero: .. code-block:: xml 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Finally, we need to assemble particle, wavefunction, and pseudopotential information into the main QMCPACK input file (``O.q0.opt.in.xml``) and specify inputs for the Jastrow optimization process. Open ``O.q0.opt.in.xml`` and write in the location of the particle, wavefunction, and pseudopotential files ("````" are comments): .. code-block:: xml ... ... ... ... The relevant portion of the input describing the linear optimization process is .. code-block:: xml ECOST UVCOST RVCOST TS SAMPLES 50 200 1 ... 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. energy Fraction of trial energy in the cost function. unreweightedvariance Fraction of unreweighted trial variance in the cost function. Neglecting the weights can be more robust. reweightedvariance Fraction of trial variance (including the full weights) in the cost function. 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 :math:`\textrm{Ha}^{-1}`. 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. 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. 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 ``qmca``). 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. 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. 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 ``loop`` section can be duplicated with a different cost function each time). Replace ``MAX``, ``EVCOST``, ``UVCOST``, ``RVCOST``, ``TS``, and ``SAMPLES`` in the ``loop`` with appropriate starting values in the ``O.q0.opt.in.xml`` input file. Perform the optimization run by typing :: mpirun -np 4 qmcpack O.q0.opt.in.xml >&O.q0.opt.out& .. jobrun_vesta qmcpack O.q0.opt.in.xml The run should take only a few minutes for reasonable values of loop ``max`` and ``samples``. The log file output will appear in ``O.q0.opt.out``. The beginning of each linear optimization will be marked with text similar to :: ========================================================= Start QMCFixedSampleLinearOptimize File Root O.q0.opt.s011 append = no ========================================================= 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: :: OldCost: 7.0598901869e-01 NewCost: 7.0592576381e-01 Delta Cost:-6.3254886314e-05 ... uu_0 6.9392504232e-01 1 1 ON 0 uu_1 4.9690781460e-01 1 1 ON 1 uu_2 4.0934542375e-01 1 1 ON 2 uu_3 3.7875640157e-01 1 1 ON 3 uu_4 3.7308380014e-01 1 1 ON 4 uu_5 3.5419786809e-01 1 1 ON 5 uu_6 4.3139019377e-01 1 1 ON 6 uu_7 1.9344371667e-01 1 1 ON 7 ud_0 3.9219009713e-01 1 1 ON 8 ud_1 1.2352664647e-01 1 1 ON 9 ud_2 4.4048945133e-02 1 1 ON 10 ud_3 2.1415676741e-02 1 1 ON 11 ud_4 1.5201803731e-02 1 1 ON 12 ud_5 2.3708169445e-02 1 1 ON 13 ud_6 3.4279064930e-02 1 1 ON 14 ud_7 4.3334583596e-02 1 1 ON 15 eO_0 -7.8490123937e-01 1 1 ON 16 eO_1 -6.6726618338e-01 1 1 ON 17 eO_2 -4.8753453838e-01 1 1 ON 18 eO_3 -3.0913993774e-01 1 1 ON 19 eO_4 -1.7901872177e-01 1 1 ON 20 eO_5 -8.6199000697e-02 1 1 ON 21 eO_6 -4.0601160841e-02 1 1 ON 22 eO_7 -4.1358075061e-03 1 1 ON 23 ... QMC Execution time = 2.8218972974e+01 secs The cost function should decrease during each linear optimization (``Delta cost < 0``). Try "``grep OldCost *opt.out.``" You should see something like this: :: 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 OldCost: 8.0649629434e-01 NewCost: 8.0551871147e-01 Delta Cost:-9.7758287036e-04 OldCost: 6.6821241388e-01 NewCost: 6.6797703487e-01 Delta Cost:-2.3537901148e-04 OldCost: 7.0106275099e-01 NewCost: 7.0078055426e-01 Delta Cost:-2.8219672877e-04 OldCost: 6.9538522411e-01 NewCost: 6.9419186712e-01 Delta Cost:-1.1933569922e-03 OldCost: 6.7709626744e-01 NewCost: 6.7501251165e-01 Delta Cost:-2.0837557922e-03 OldCost: 6.6659923822e-01 NewCost: 6.6651737755e-01 Delta Cost:-8.1860671682e-05 OldCost: 7.7828995609e-01 NewCost: 7.7735482525e-01 Delta Cost:-9.3513083900e-04 OldCost: 7.2717974404e-01 NewCost: 7.2715201115e-01 Delta Cost:-2.7732880747e-05 OldCost: 6.9400639873e-01 NewCost: 6.9257183689e-01 Delta Cost:-1.4345618444e-03 OldCost: 7.0598901869e-01 NewCost: 7.0592576381e-01 Delta Cost:-6.3254886314e-05 Blocked averages of energy data, including the kinetic energy and components of the potential energy, are written to ``scalar.dat`` files. The first is named "``O.q0.opt.s000.scalar.dat``," with a series number of zero (``s000``). In the end there will be ``MAX`` of them, one for each series. When the job has finished, use the ``qmca`` tool to assess the effectiveness of the optimization process. To look at just the total energy and the variance, type "``qmca -q ev O.q0.opt*scalar*``." This will print the energy, variance, and the variance/energy ratio in Hartree units: :: LocalEnergy Variance ratio O.q0.opt series 0 -15.739585 +/- 0.007656 0.887412 +/- 0.010728 0.0564 O.q0.opt series 1 -15.848347 +/- 0.004089 0.318490 +/- 0.006404 0.0201 O.q0.opt series 2 -15.867494 +/- 0.004831 0.292309 +/- 0.007786 0.0184 O.q0.opt series 3 -15.871508 +/- 0.003025 0.275364 +/- 0.006045 0.0173 O.q0.opt series 4 -15.865512 +/- 0.002997 0.278056 +/- 0.006523 0.0175 O.q0.opt series 5 -15.864967 +/- 0.002733 0.278065 +/- 0.004413 0.0175 O.q0.opt series 6 -15.869644 +/- 0.002949 0.273497 +/- 0.006141 0.0172 O.q0.opt series 7 -15.868397 +/- 0.003838 0.285451 +/- 0.007570 0.0180 ... Plots of the data can also be obtained with the “``-p``” option (“``qmca -p -q ev O.q0.opt*scalar*``”). Identify which optimization series is the "best" according to your cost function. It is likely that multiple series are similar in quality. Note the ``opt.xml`` file corresponding to this series. This file contains the final value of the optimized Jastrow parameters to be used in the DMC calculations of the next section of the lab. **Questions and Exercises** #. What is the acceptance ratio of your optimization runs? (use “texttqmca -q ar O.q0.opt*scalar\*”) Do you expect the MC sampling to be efficient? #. How do you know when the optimization process has converged? #. (optional) Optimization is sometimes sensitive to initial guesses of the parameters. If you have time, try varying the initial parameters, including the cutoff radius (``rcut``) of the Jastrow factors (remember to change ``id`` in the ```` element). Do you arrive at a similar set of final Jastrow parameters? What is the lowest variance you are able to achieve? 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. In the same directory you used to perform wavefunction optimization (``oxygen_atom``) you will find a sample DMC input file for the neutral oxygen atom named ``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 ``pw2qmcpack`` but instead should come from the optimization run: .. code-block:: xml Replace "``OPT_XML``" with the ``opt.xml`` file corresponding to the best Jastrow parameters you found in the last section (this is a file name similar to ``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. Descriptions of the input keywords follow. 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%. warmupSteps Number of MC steps discarded as a warmup or equilibration period of the random walk. steps Number of MC steps per block. Physical quantities, such as the total energy, are averaged over walkers and steps. 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 ``qmca``). The total number of MC steps each walker takes is ``blocks``\ :math:`\times`\ ``steps``. 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. 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). The purpose of the VMC run is to provide initial electron positions for each DMC walker. Setting :math:`\texttt{walkers}=1` in the VMC block ensures there will be only one VMC walker per execution thread. There will be a total of 4 VMC walkers in this case (see ``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 :math:`^{-1}` in imaginary time (in general, this is system dependent). This leads to a rough rule of thumb for choosing ``blocks`` and ``steps`` for the VMC run (:math:`\texttt{vwalkers}=4` here): .. math:: :label: eq68 \begin{aligned} \texttt{VBLOCKS}\times\texttt{VSTEPS} \ge \frac{\texttt{DWALKERS}}{\texttt{VWALKERS}} \frac{5~\textrm{Ha}^{-1}}{\texttt{VTIMESTEP}}\end{aligned} 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 :math:`\sim 2,000` walkers, although the dependence of the total energy on population size should be explicitly checked from time to time. To study time step bias, we will perform a sequence of DMC runs over a range of time steps (:math:`0.1` Ha\ :math:`^{-1}` is too large, and time steps below :math:`0.002` Ha\ :math:`^{-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 .. math:: :label: eq69 \begin{aligned} \texttt{timestep}_{n} &= \texttt{timestep}_{n-1}/2\nonumber\\ \texttt{warmupSteps}_{n} &= \texttt{warmupSteps}_{n-1}\times 2\nonumber\\ \texttt{blocks}_{n} &= \texttt{blocks}_{n-1}\nonumber\\ \texttt{steps}_{n} &= \texttt{steps}_{n-1}\times 2\end{aligned} 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. :math:`\texttt{blocks}\times\texttt{steps}\times\texttt{timestep}\sim 60~\mathrm{Ha}^{-1}` is sufficient for this system. Choose an initial DMC time step and create a sequence of :math:`N` time steps according to :eq:`eq69`. Make :math:`N` copies of the DMC XML block in the input file. .. code-block:: xml DWARMUP DBLOCKS DSTEPS DTIMESTEP yes Fill in ``DWARMUP``, ``DBLOCKS``, ``DSTEPS``, and ``DTIMESTEP`` for each DMC run according to :eq:`eq69`. Start the DMC time step extrapolation run by typing: :: mpirun -np 4 qmcpack O.q0.dmc.in.xml >&O.q0.dmc.out& .. else jobrun_vesta qmcpack O.q0.dmc.in.xml The run should take only a few minutes to complete. QMCPACK will create files prefixed with ``O.q0.dmc``. The log file is ``O.q0.dmc.out``. As before, block-averaged data is written to ``scalar.dat`` files. In addition, DMC runs produce ``dmc.dat`` files, which contain energy data averaged only over the walker population (one line per DMC step). The ``dmc.dat`` files also provide a record of the walker population at each step. .. _fig16: .. figure:: /figs/lab_qmc_basics_timestep_conv.jpg :width: 500 :align: center Linear fit to DMC timestep data from ``PlotTstepConv.pl``. Use the ``PlotTstepConv.pl`` to obtain a linear fit to the time step data (type "``PlotTstepConv.pl O.q0.dmc.in.xml 40``"). You should see a plot similar to :numref:`fig16`. The tail end of the text output displays the parameters for the linear fit. The "``a``" parameter is the total energy extrapolated to zero time step in Hartree units. :: ... Final set of parameters Asymptotic Standard Error ======================= ========================== a = -15.8925 +/- 0.0007442 (0.004683%) b = -0.0457479 +/- 0.0422 (92.24%) ... **Questions and Exercises** #. What is the :math:`\tau\rightarrow 0` extrapolated value for the total energy? #. What is the maximum time step you should use if you want to calculate the total energy to an accuracy of :math:`0.05` eV? For convenience, :math:`1~\textrm{Ha}=27.2113846~\textrm{eV}`. #. What is the acceptance ratio for this (bias :math:`<0.05` eV) run? Does it follow the rule of thumb for sensible DMC (acceptance ratio :math:`>99`\ %) ? #. Check the fluctuations in the walker population (``qmca -t -q nw O.q0.dmc*dmc.dat –noac``). Does the population seem to be stable? #. (Optional) Study population control bias for the oxygen atom. Select a few population sizes. Copy ``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 “textttsamples,” and choose a unique ``id`` in ````. Use ``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? 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 :math:`+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 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. 1. Obtain DFT orbitals with QE. (a) Copy the DFT input (``O.q0.dft.in``) to ``O.q1.dft.in`` (b) Edit ``O.q1.dft.in`` to match the +1 charge state of the oxygen atom. :: ... prefix = 'O.q1' ... tot_charge = 1 tot_magnetization = 3 ... (c) Perform the DFT run: ``mpirun -np 4 pw.x -input O.q1.dft.in >&O.q1.dft.out&`` 2. Convert the orbitals to ESHDF format. (a) Copy the pw2qmcpack input (``O.q0.p2q.in``) to ``O.q1.p2q.in`` (b) Edit ``O.q1.p2q.in`` to match the file prefix used in DFT. :: ... prefix = 'O.q1' ... (c) Perform the orbital conversion run: ``mpirun -np 1 pw2qmcpack.x&O.q1.p2q.out&`` 3. Optimize the Jastrow factor with QMCPACK. (a) Copy the optimization input (``O.q0.opt.in.xml``) to ``O.q1.opt.in.xml`` (b) Edit ``O.q1.opt.in.xml`` to match the file prefix used in DFT. .. code-block:: xml ... ... ... ... (c) Edit the particle XML file (``O.q1.ptcl.xml``) to have open boundary conditions. .. code-block:: xml n n n (d) Add cutoffs to the Jastrow factors in the wavefunction XML file (``O.q1.wfs.xml``) .. code-block:: xml ... ... ... ... (e) Perform the Jastrow optimization run: ``mpirun -np 4 qmcpack O.q1.opt.in.xml >&O.q1.opt.out&`` (f) Identify the optimal set of parameters with ``qmca`` (``[your opt.xml]``). 4. DMC time step study with QMCPACK (a) Copy the DMC input (``O.q0.dmc.in.xml``) to ``O.q1.dmc.in.xml`` (b) Edit ``O.q1.dmc.in.xml`` to use the DFT prefix and the optimal Jastrow. .. code-block:: xml ... ... ... ... (c) Perform the DMC run: ``mpirun -np 4 qmcpack O.q1.dmc.in.xml >&O.q1.dmc.out&`` (d) Obtain the DMC total energy extrapolated to zero time step with ``PlotTstepConv.pl``. 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 :ref:`dimer-automation` of the lab. **Questions and Exercises** #. What is the :math:`\tau\rightarrow 0` extrapolated DMC value for the first ionization potential of oxygen? #. How does the extrapolated value compare with the experimental IP? Go to http://physics.nist.gov/PhysRefData/ASD/ionEnergy.html and enter “``O I``” in the box labeled “``Spectra``” and click on the “``Retrieve Data``” button. #. What can we conclude about the accuracy of the pseudopotential? What factors complicate this assessment? #. Explore the sensitivity of the IP to the choice of time step. Type ``./ip_conv.py`` to view three time step extrapolation plots: two for the :math:`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? #. What is the maximum time step you should use if you want to calculate the ionization potential to an accuracy of :math:`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? #. Are the acceptance ratio and population fluctuations reasonable for the :math:`q=1` calculations? 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: #. Orbital generation via QE or GAMESS. #. Conversion of orbital data via ``pw2qmcpack.x`` or ``convert4qmc``. #. Optimization of Jastrow factors via QMCPACK. #. DMC calculation via QMCPACK. 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 (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 pregenerated orbitals) follows. Take a moment to read it over and especially note the comments (prefixed with "``\#``") explaining most of the contents. If the input syntax is unclear you may want to consult portions of :ref:`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`. .. code-block:: python #! /usr/bin/env python3 # import Nexus functions from nexus import settings,job,get_machine,run_project from nexus import generate_physical_system from nexus import generate_qmcpack,vmc settings( # Nexus settings pseudo_dir = './pseudopotentials', # location of PP files runs = '', # root directory for simulations results = '', # root directory for simulation results status_only = 0, # show simulation status, then exit generate_only = 0, # generate input files, then exit sleep = 3, # seconds between checks on sim. progress machine = 'ws4', # workstation with 4 cores ) qmcjob = job( # specify job parameters cores = 4, # use 4 MPI tasks threads = 1, # 1 OpenMP thread per node app = 'qmcpack' # use QMCPACK executable (assumed in PATH) ) qmc_calcs = [ # list QMC calculation methods vmc( # VMC walkers = 1, # 1 walker warmupsteps = 50, # 50 MC steps for warmup blocks = 200, # 200 blocks steps = 10, # 10 steps per block timestep = .4 # 0.4 1/Ha timestep )] dimer = generate_physical_system( # make a dimer system type = 'dimer', # system type is dimer dimer = ('O','O'), # dimer is two oxygen atoms separation = 1.2074, # separated by 1.2074 Angstrom Lbox = 15.0, # simulation box is 15 Angstrom units = 'A', # Angstrom is dist. unit net_spin = 2, # nup-ndown is 2 O = 6 # pseudo-oxygen has 6 valence el. ) qmc = generate_qmcpack( # make a qmcpack simulation identifier = 'example', # prefix files with 'example' path = 'scale_1.0', # run in ./scale_1.0 directory system = dimer, # run the dimer system job = qmcjob, # set job parameters input_type = 'basic', # basic qmcpack inputs given below pseudos = ['O.BFD.xml'], # list of PP's to use orbitals_h5 = 'O2.pwscf.h5', # file with orbitals from DFT bconds = 'nnn', # open boundary conditions jastrows = [], # no jastrow factors calculations = qmc_calcs # QMC calculations to perform ) run_project(qmc) # write input file and submit job .. _dimer-automation: Automated binding curve of the oxygen dimer ------------------------------------------- In this section we will use Nexus to calculate the DMC total energy of the oxygen dimer over a series of bond lengths. The equilibrium bond length and binding energy of the dimer will be determined by performing a polynomial fit to the data (Morse potential fits should be preferred in production tests). Comparing these values with corresponding experimental data provides a second test of the BFD pseudopotential for oxygen. Enter the ``oxygen_dimer`` directory. Copy your BFD pseudopotential from the atom runs into ``oxygen_dimer/pseudopotentials`` (be sure to move both files: ``.upf`` and ``.xml``). Open ``O_dimer.py`` with a text editor. The overall format is similar to the example file shown in the last section. The main difference is that a full workflow of runs (DFT orbital generation, orbital conversion, optimization and DMC) are being performed rather than a single VMC run. As in the example in the last section, the oxygen dimer is generated with the ``generate_physical_ system`` function: .. code-block:: python dimer = generate_physical_system( type = 'dimer', dimer = ('O','O'), separation = 1.2074*scale, Lbox = 10.0, units = 'A', net_spin = 2, O = 6 ) Similar syntax can be used to generate crystal structures or to specify systems with arbitrary atomic configurations and simulation cells. Notice that a "``scale``" variable has been introduced to stretch or compress the dimer. Next, objects representing a QE (PWSCF) run and subsequent orbital conversion step are constructed with respective ``generate_*`` functions: .. code-block:: python dft = generate_pwscf( identifier = 'dft', ... input_dft = 'lda', ... ) sims.append(dft) # describe orbital conversion run p2q = generate_pw2qmcpack( identifier = 'p2q', ... dependencies = (dft,'orbitals'), ) sims.append(p2q) Note the ``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 ``generate_qmcpack`` function: .. code-block:: python opt = generate_qmcpack( identifier = 'opt', ... jastrows = [('J1','bspline',8,5.0), ('J2','bspline',8,10.0)], calculations = [ loop(max=12, qmc=linear( energy = 0.0, unreweightedvariance = 1.0, reweightedvariance = 0.0, timestep = 0.3, samples = 61440, warmupsteps = 50, blocks = 200, substeps = 1, walkers = 1, minwalkers = 0.5, maxweight = 1e9, usedrift = False, minmethod = 'quartic', beta = 0.025, exp0 = -16, bigchange = 15.0, alloweddifference = 1e-4, stepsize = 0.2, stabilizerscale = 1.0, nstabilizers = 3, ) ) ], dependencies = (p2q,'orbitals'), ) sims.append(opt) qmc = generate_qmcpack( identifier = 'qmc', ... jastrows = [], calculations = [ vmc( walkers = 1, warmupsteps = 30, blocks = 20, steps = 10, substeps = 2, timestep = .4, samples = 2048 ), dmc( warmupsteps = 100, blocks = 400, steps = 32, timestep = 0.01, nonlocalmoves = True, ) ], dependencies = [(p2q,'orbitals'),(opt,'jastrow')], ) sims.append(qmc) Shared details such as the run directory, job, pseudopotentials, and orbital file have been omitted (``...``). The "``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 ``dependencies`` indicates that the DMC run depends on the optimization run for the Jastrow factor. Nexus will submit the "``opt``" run first, and upon completion it will scan the output, select the optimal set of parameters, pass the Jastrow information to the "``qmc``" run, and then submit the DMC job. Independent job workflows are submitted in parallel when permitted. No input files are written or job submissions made until the "``run_project``" function is reached: .. code-block:: python run_project(sims) All of the simulation objects have been collected into a list (``sims``) for submission. As written, ``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 (``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 (``scale<1.0``) and stretch (``scale>1.0``): .. code-block:: python scales = [1.00,0.90,0.95,1.05,1.10] Note that "``1.00``" is left in front because we are going to optimize the Jastrow factor first at the equilibrium separation and reuse this Jastrow factor for all other separation distances. This procedure is used because it can reduce variations in localization errors (due to pseudopotentials in DMC) along the binding curve. Change the ``status_only`` parameter in the "``settings``" function to ``1`` and type "``./O_dimer.py``" at the command line. This will print the status of all simulations: :: Project starting checking for file collisions loading cascade images cascade 0 checking in cascade 10 checking in cascade 4 checking in cascade 13 checking in cascade 7 checking in checking cascade dependencies all simulation dependencies satisfied cascade status setup, sent_files, submitted, finished, got_output, analyzed 000000 dft ./scale_1.0 000000 p2q ./scale_1.0 000000 opt ./scale_1.0 000000 qmc ./scale_1.0 000000 dft ./scale_0.9 000000 p2q ./scale_0.9 000000 qmc ./scale_0.9 000000 dft ./scale_0.95 000000 p2q ./scale_0.95 000000 qmc ./scale_0.95 000000 dft ./scale_1.05 000000 p2q ./scale_1.05 000000 qmc ./scale_1.05 000000 dft ./scale_1.1 000000 p2q ./scale_1.1 000000 qmc ./scale_1.1 setup, sent_files, submitted, finished, got_output, analyzed In this case, five simulation "cascades" (workflows) have been identified, each one starting and ending with "``dft``" and "``qmc``" runs, respectively. The six status flags ``setup``, ``sent_files``, ``submitted``, ``finished``, ``got_output``, ``analyzed``) each shows ``0``, indicating that no work has been done yet. Now change "``status_only``" back to ``0``, set "``generate_only``" to ``1``, and run ``O_dimer.py`` again. This will perform a dry run of all simulations. The dry run should finish in about 20 seconds: :: Project starting checking for file collisions loading cascade images cascade 0 checking in cascade 10 checking in cascade 4 checking in cascade 13 checking in cascade 7 checking in checking cascade dependencies all simulation dependencies satisfied starting runs: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ poll 0 memory 91.03 MB Entering ./scale_1.0 0 writing input files 0 dft Entering ./scale_1.0 0 sending required files 0 dft submitting job 0 dft ... poll 1 memory 91.10 MB ... Entering ./scale_1.0 0 Would have executed: export OMP_NUM_THREADS=1 mpirun -np 4 pw.x -input dft.in poll 2 memory 91.10 MB Entering ./scale_1.0 0 copying results 0 dft Entering ./scale_1.0 0 analyzing 0 dft ... poll 3 memory 91.10 MB Entering ./scale_1.0 1 writing input files 1 p2q Entering ./scale_1.0 1 sending required files 1 p2q submitting job 1 p2q ... Entering ./scale_1.0 1 Would have executed: export OMP_NUM_THREADS=1 mpirun -np 1 pw2qmcpack.x&scf.out& (wait until JOB DONE appears in scf.out) mpirun -np 1 pw2qmcpack.x&p2q.out& Once the conversion process has finished, the orbitals should be located in the file ``diamond_vmc/pwscf_output/pwscf.pwscf.h5``. Open ``diamond_vmc/vmc.in.xml`` and replace "``MISSING.h5``" with "``./pwscf_output/pwscf.pwscf.h5``". Next submit the VMC run: :: mpirun -np 4 qmcpack vmc.in.xml>&vmc.out& 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 ``example.py`` in the commented out text. .. _python-basics: Appendix A: Basic Python constructs ----------------------------------- Basic Python data types (``int``, ``float``, ``str``, ``tuple``, ``list``, ``array``, ``dict``, ``obj``) and programming constructs (``if`` statements, ``for`` loops, functions w/ keyword arguments) are briefly overviewed in the following. All examples can be executed interactively in Python. To do this, type "``python``" at the command line and paste any of the shaded text below at the ``>>>`` prompt. For more information about effective use of Python, consult the detailed online documentation: https://docs.python.org/2/. Intrinsic types: ``int, float, str`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: python #this is a comment i=5 # integer f=3.6 # float s='quantum/monte/carlo' # string n=None # represents "nothing" f+=1.4 # add-assign (-,*,/ also): 5.0 2**3 # raise to a power: 8 str(i) # int to string: '5' s+'/simulations' # joining strings: 'quantum/monte/carlo/simulations' 'i={0}'.format(i) # format string: 'i=5' Container types: ``tuple, list, array, dict, obj`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: python from numpy import array # get array from numpy module from developer import obj # get obj from Nexus' developer module t=('A',42,56,123.0) # tuple l=['B',3.14,196] # list a=array([1,2,3]) # array d={'a':5,'b':6} # dict o=obj(a=5,b=6) # obj # printing print(t) # ('A', 42, 56, 123.0) print(l) # ['B', 3.1400000000000001, 196] print(a) # [1 2 3] print(d) # {'a': 5, 'b': 6} print(o) # a = 5 # b = 6 len(t),len(l),len(a),len(d),len(o) #number of elements: (4, 3, 3, 2, 2) t[0],l[0],a[0],d['a'],o.a #element access: ('A', 'B', 1, 5, 5) s = array([0,1,2,3,4]) # slices: works for tuple, list, array s[:] # array([0, 1, 2, 3, 4]) s[2:] # array([2, 3, 4]) s[:2] # array([0, 1]) s[1:4] # array([1, 2, 3]) s[0:5:2] # array([0, 2, 4]) # list operations l2 = list(l) # make independent copy l.append(4) # add new element: ['B', 3.14, 196, 4] l+[5,6,7] # addition: ['B', 3.14, 196, 4, 5, 6, 7] 3*[0,1] # multiplication: [0, 1, 0, 1, 0, 1] b=array([5,6,7]) # array operations a2 = a.copy() # make independent copy a+b # addition: array([ 6, 8, 10]) a+3 # addition: array([ 4, 5, 6]) a*b # multiplication: array([ 5, 12, 21]) 3*a # multiplication: array([3, 6, 9]) # dict/obj operations d2 = d.copy() # make independent copy d['c'] = 7 # add/assign element d.keys() # get element names: ['a', 'c', 'b'] d.values() # get element values: [5, 7, 6] # obj-specific operations o.c = 7 # add/assign element o.set(c=7,d=8) # add/assign multiple elements 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 ``obj`` instance in the following example, but it also holds for ``list``, ``array``, ``dict``, and others. .. code-block:: python >>> o = obj(a=5,b=6) >>> >>> p=o >>> >>> p.a=7 >>> >>> print(o) a = 7 b = 6 >>> q=o.copy() >>> >>> q.a=9 >>> >>> print(o) a = 7 b = 6 Here ``p`` is just another name for ``o``, while ``q`` is a fully independent copy of it. Conditional Statements: ``if/elif/else`` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: python a = 5 if a is None: print('a is None') elif a==4: print('a is 4') elif a<=6 and a>2: print('a is in the range (2,6]') elif a<-1 or a>26: print('a is not in the range [-1,26]') elif a!=10: print('a is not 10') else: print('a is 10') #end if The "``\#end if``" is not part of Python syntax, but you will see text like this throughout Nexus for clear encapsulation. Iteration: ``for`` ^^^^^^^^^^^^^^^^^^ .. code-block:: python from developer import obj l = [1,2,3] m = [4,5,6] s = 0 for i in range(len(l)): # loop over list indices s += l[i] + m[i] #end for print(s) # s is 21 s = 0 for v in l: # loop over list elements s += v #end for print(s) # s is 6 o = obj(a=5,b=6) s = 0 for v in o: # loop over obj elements s += v #end for print(s) # s is 11 d = {'a':5,'b':4} for n,v in o.items():# loop over name/value pairs in obj d[n] += v #end for print(d) # d is {'a': 10, 'b': 10} Functions: ``def``, argument syntax ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code-block:: python def f(a,b,c=5): # basic function, c has a default value print(a,b,c) #end def f f(1,b=2) # prints: 1 2 5 def f(*args,**kwargs): # general function, returns nothing print(args) # args: tuple of positional arguments print(kwargs) # kwargs: dict of keyword arguments #end def f f('s',(1,2),a=3,b='t') # 2 pos., 2 kw. args, prints: # ('s', (1, 2)) # {'a': 3, 'b': 't'} l = [0,1,2] f(*l,a=6) # pos. args from list, 1 kw. arg, prints: # (0, 1, 2) # {'a': 6} o = obj(a=5,b=6) f(*l,**o) # pos./kw. args from list/obj, prints: # (0, 1, 2) # {'a': 5, 'b': 6} f( # indented kw. args, prints blocks = 200, # () steps = 10, # {'steps': 10, 'blocks': 200, 'timestep': 0.01} timestep = 0.01 ) o = obj( # obj w/ indented kw. args blocks = 100, steps = 5, timestep = 0.02 ) f(**o) # kw. args from obj, prints: # () # {'timestep': 0.02, 'blocks': 100, 'steps': 5} .. bibliography:: /bibs/labs_qmc_basics.bib