Periodic LCAO for Solids¶
Introduction¶
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
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.
#! /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
cell.build()
sp_twist=[0.07761248, 0.07761248, -0.07761248]
kmesh=[2,1,1]
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
e_scf=mf.kernel()
ener = open('e_scf','w')
ener.write('%s\n' % (e_scf))
print('e_scf',e_scf)
ener.close()
title="C_diamond-tiled-cplx"
from PyscfToQmcpack import savetoqmcpack
savetoqmcpack(cell,mf,title=title,kmesh=kmesh,kpts=kpts,sp_twist=sp_twist)
Note that the last three lines of the file
title="C_diamond-tiled-cplx"
from PyscfToQmcpack import savetoqmcpack
savetoqmcpack(cell,mf,title=title,kmesh=kmesh,kpts=kpts,sp_twist=sp_twist)
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 PyscfToQmcpack.py. 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: https://github.com/QMCPACK/qmcpack_workshop_2019 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
export PYTHONPATH=QMCPACK_PATH/src/QMCTools:$PYTHONPATH
To generate QMCPACK input files, you will need to run convert4qmc
exactly as specified in convert4qmc for both cases:
convert4qmc -pyscf 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:
<?xml version="1.0"?>
<qmcsystem>
<simulationcell>
<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>
<parameter name="bconds">p p p</parameter>
<parameter name="LR_dim_cutoff">15</parameter>
</simulationcell>
<particleset name="ion0" size="4">
<group name="C">
<parameter name="charge">4</parameter>
<parameter name="valence">4</parameter>
<parameter name="atomicnumber">6</parameter>
</group>
<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>
<attrib name="ionid" datatype="stringArray">
C C C C
</attrib>
</particleset>
<particleset name="e" random="yes" randomsrc="ion0">
<group name="u" size="8">
<parameter name="charge">-1</parameter>
</group>
<group name="d" size="8">
<parameter name="charge">-1</parameter>
</group>
</particleset>
</qmcsystem>
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.
<?xml version="1.0"?>
<qmcsystem>
<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">
<slaterdeterminant>
<determinant id="updet" size="8">
<occupation mode="ground"/>
<coefficient size="52" spindataset="0"/>
</determinant>
<determinant id="downdet" size="8">
<occupation mode="ground"/>
<coefficient size="52" spindataset="0"/>
</determinant>
</slaterdeterminant>
</determinantset>
<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>
<correlation size="10" speciesA="u" speciesB="d">
<coefficients id="ud" type="Array"> 0 0 0 0 0 0 0 0 0 0</coefficients>
</correlation>
</jastrow>
<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>
</correlation>
</jastrow>
</wavefunction>
</qmcsystem>
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
tag type
default
description
twist
3 doubles
(0 0 0)
Coordinate of the twist to compute
href
string
default
Name of the HDF5 file generated by PySCF and used for convert4qmc
PBCimages
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.
- SBB+18(1,2)
Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D. McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan. Pyscf: the python-based simulations of chemistry framework. Wiley Interdisciplinary Reviews: Computational Molecular Science, 8(1):n/a–n/a, 2018. doi:10.1002/wcms.1340.