Periodic LCAO for Solids


QMCPACK implements the linear combination of atomic orbitals (LCAO) and Gaussian basis sets in periodic boundary conditions. This method uses orders of magnitude less memory than the real-space spline wavefunction. Although the spline scheme enables very fast evaluation of the wavefunction, it might require too much on-node memory for a large complex cell. The periodic Gaussian evaluation provides a fallback that will definitely fit in available memory but at significantly increased computational expense. Well-designed Gaussian basis sets should be used to accurately represent the wavefunction, typically including both diffuse and high angular momentum functions.

The current implementation is not highly optimized for efficiency but can handle real and complex trial wavefunctions generated by PySCF [SBB+18], but other codes such as Crystal can be interfaced on request. Supercell tiling is handled outside QMCPACK through a proper PySCF input generated by Nexus and the Supercell geometry and coefficients of the molecular orbotals are constructed in the converter provided by QMCPACK. This is different from the plane wave/spline route where the tiling is provided in QMCPACK.

LCAO schemes use physical considerations to construct a highly efficient basis set compared with plane waves. Typically only a few tens of basis functions per atom are required compared with thousands of plane waves. Many forms of LCAO schemes exist and are being implemented in QMCPACK. The details of the already-implemented methods are described in the following section.

GTOs: The Gaussian basis functions follow a radial-angular decomposition of

(56)\[ \phi (\mathbf{r} )=R_{l}(r)Y_{lm}(\theta ,\phi )\:,\]

where \(Y_{{lm}}(\theta ,\phi )\) is a spherical harmonic, \(l\) and \(m\) are the angular momentum and its \(z\) component, and \(r, \theta, \phi\) are spherical coordinates. In practice, they are atom centered and the \(l\) expansion typically includes 1–3 additional channels compared with the formally occupied states of the atom (e.g., 4–6 for a nickel atom with occupied \(s\), \(p\), and \(d\) electron shells.

The evaluation of GTOs within PBC differs slightly from evaluating GTOs in open boundary conditions (OBCs). The orbitals are evaluated at a distance \(r\) in the primitive cell (similar to OBC), and then the contributions of the periodic images are added by evaluating the orbital at a distance \(r+T\), where T is a translation of the cell lattice vector. This requires loops over the periodic images until the contributions are orbitals \(\Phi\). In the current implementation, the number of periodic images is an input parameter named PBCimages, which takes three integers corresponding to the number of periodic images along the supercell axes (X, Y and Z axes for a cubic cell). By default these parameters are set to PBCimages= 8 8 8, but they require manual convergence checks. Convergence checks can be performed by checking the total energy convergence with respect to PBCimages, similar to checks performed for plane wave cutoff energy and B-spline grids. Use of diffuse Gaussians might require these parameters to be increased, while sharply localized Gaussians might permit a decrease. The cost of evaluating the wavefunction increases sharply as PBCimages is increased. This input parameter will be replaced by a tolerance factor and numerical screening in the future.

Generating and using periodic Gaussian-type wavefunctions using PySCF

Similar to any QMC calculation, using periodic GTOs requires the generation of a periodic trial wavefunction. QMCPACK is currently interfaced to PySCF, which is a multipurpose electronic structure written mainly in Python with key numerical functionality implemented via optimized C and C++ libraries [SBB+18]. Such a wavefunction can be generated according to the following example for a \(2 \times 1 \times 1\) supercell using tiling (kpoints) and a supertwist shifted away from \(\Gamma\), leading to a complex wavefunction.

Listing 57 Example PySCF input for single k-point calculation for a \(2 \times 1 \times 1\) carbon supercell.
#! /usr/bin/env python3
import numpy
import h5py
from pyscf.pbc import gto, scf, dft, df
from pyscf.pbc import df

cell = gto.Cell()
cell.a             = '''
         3.37316115       3.37316115       0.00000000
         0.00000000       3.37316115       3.37316115
         3.37316115       0.00000000       3.37316115'''
cell.atom = '''
   C        0.00000000       0.00000000       0.00000000
   C        1.686580575      1.686580575      1.686580575
cell.basis         = 'bfd-vdz'
cell.ecp           = 'bfd'
cell.unit          = 'B'
cell.drop_exponent = 0.1
cell.verbose       = 5
cell.charge        = 0
cell.spin          = 0

sp_twist=[0.07761248, 0.07761248, -0.07761248]

kpts=[[ 0.07761248,  0.07761248, -0.07761248],[ 0.54328733,  0.54328733, -0.54328733]]

mf = scf.KRHF(cell,kpts)
mf.exxdiv = 'ewald'
mf.max_cycle = 200


ener = open('e_scf','w')
ener.write('%s\n' % (e_scf))

from PyscfToQmcpack import savetoqmcpack

Note that the last three lines of the file

from PyscfToQmcpack import savetoqmcpack

contain the title (name of the HDF5 to be used in QMCPACK) and the call to the converter. The title variable will be the name of the HDF5 file where all the data needed by QMCPACK will be stored. The function savetoqmcpack will be called at the end of the calculation and will generate the HDF5 similarly to the nonperiodic PySCF calculation in convert4qmc. The function is distributed with QMCPACK and is located in the qmcpack/src/QMCTools directory under the name Note that you need to specify the supertwist coordinates that was used with the provided kpoints. The supertwist must match the coordinates of the K-points otherwise the phase factor for the atomic orbital will be incorrect and incorrect results will be obtained. (For more details on how to generate tiling with PySCF and Nexus, refer to the Nexus guide or the 2019 QMCPACK Workshop material available on github: under qmcpack_workshop_2019/day2_nexus/pyscf/04_pyscf_diamond_hf_qmc/

For the converter in the script to be called properly, you need to specify the path to the file in your PYTHONPATH such as


To generate QMCPACK input files, you will need to run convert4qmc exactly as specified in convert4qmc for both cases:

convert4qmc -orbitals C_diamond-tiled-cplx

This tool can be used with any option described in convert4qmc. Since the HDF5 contains all the information needed, there is no need to specify any other specific tag for periodicity. A supercell at \(\Gamma\)-point or using multiple k-points will work without further modification.

Running convert4qmc will generate 3 input files:

Listing 58 C_diamond-tiled-cplx.structure.xml. This file contains the geometry of the system.
<?xml version="1.0"?>
    <parameter name="lattice">
  6.74632230000000e+00  6.74632230000000e+00  0.00000000000000e+00
  0.00000000000000e+00  3.37316115000000e+00  3.37316115000000e+00
  3.37316115000000e+00  0.00000000000000e+00  3.37316115000000e+00
    <parameter name="bconds">p p p</parameter>
    <parameter name="LR_dim_cutoff">15</parameter>
  <particleset name="ion0" size="4">
    <group name="C">
      <parameter name="charge">4</parameter>
      <parameter name="valence">4</parameter>
      <parameter name="atomicnumber">6</parameter>
    <attrib name="position" datatype="posArray">
  0.0000000000e+00  0.0000000000e+00  0.0000000000e+00
  1.6865805750e+00  1.6865805750e+00  1.6865805750e+00
  3.3731611500e+00  3.3731611500e+00  0.0000000000e+00
  5.0597417250e+00  5.0597417250e+00  1.6865805750e+00
    <attrib name="ionid" datatype="stringArray">
 C C C C
  <particleset name="e" random="yes" randomsrc="ion0">
    <group name="u" size="8">
      <parameter name="charge">-1</parameter>
    <group name="d" size="8">
      <parameter name="charge">-1</parameter>

As one can see, for both examples, the two-atom primitive cell has been expanded to contain four atoms in a \(2 \times 1 \times 1\) carbon cell.

Listing 59 C_diamond-tiled-cplx.wfj.xml. This file contains the trial wavefunction.
<?xml version="1.0"?>
  <wavefunction name="psi0" target="e">
    <determinantset type="MolecularOrbital" name="LCAOBSet" source="ion0" transform="yes" twist="0.07761248  0.07761248  -0.07761248" href="C_diamond-tiled-cplx.h5" PBCimages="8  8  8">
        <determinant id="updet" size="8">
          <occupation mode="ground"/>
          <coefficient size="52" spindataset="0"/>
        <determinant id="downdet" size="8">
          <occupation mode="ground"/>
          <coefficient size="52" spindataset="0"/>

    <jastrow name="J2" type="Two-Body" function="Bspline" print="yes">
      <correlation size="10" speciesA="u" speciesB="u">
        <coefficients id="uu" type="Array"> 0 0 0 0 0 0 0 0 0 0</coefficients>
      <correlation size="10" speciesA="u" speciesB="d">
        <coefficients id="ud" type="Array"> 0 0 0 0 0 0 0 0 0 0</coefficients>
    <jastrow name="J1" type="One-Body" function="Bspline" source="ion0" print="yes">
      <correlation size="10" cusp="0" elementType="C">
        <coefficients id="eC" type="Array"> 0 0 0 0 0 0 0 0 0 0</coefficients>

This file contains information related to the trial wavefunction. It is identical to the input file from an OBC calculation to the exception of the following tags:

*.wfj.xml specific tags:


tag type




3 doubles

(0 0 0)

Coordinate of the twist to compute




Name of the HDF5 file generated by PySCF and used for convert4qmc


3 Integer

8 8 8

Number of periodic images to evaluate the orbitals

Other files containing QMC methods (such as optimization, VMC, and DMC blocks) will be generated and will behave in a similar fashion regardless of the type of SPO in the trial wavefunction.