Benzene with a Jastrow–Slater single-determinant (JSD) and Jastrow–Antisymmetrized Geminal Function ansatz (JAGP)

00 Introduction 

In this tutorial, you will learn a complete workflow for VMC and LRDMC calculations on the benzene molecule with JSD and JAGP ansatz, starting from a DFT calculation in PySCF that employs correlation-consistent effective core potentials (ccECPs). All input and output files for this tutorial can be downloaded here.

01 PySCF calculation and its conversion to a TREXIO file 

The first step is to run a PySCF calculation by typing:

# pyscf calculation
cd 01_pyscf_calculation
python3 pyscf_benzene.py

The Python code is given as follows:

#!/usr/bin/env python
# coding: utf-8

# pySCF -> pyscf checkpoint file (NH3 molecule)

# load python packages
import os, sys

# load ASE modules
from ase.io import read

# load pyscf packages
from pyscf import gto, scf, mp, tools
 
#open boundary condition
structure_file="benzene.xyz"
checkpoint_file="benzene.chk"
pyscf_output="out_benzene_pyscf"
charge=0
spin=0
basis="ccecp-ccpvdz"
ecp='ccecp'
scf_method="HF"  # HF or DFT
dft_xc="LDA_X,LDA_C_PZ" # XC for DFT

print(f"structure file = {structure_file}")

# read a structure
atom=read(structure_file)
chemical_symbols=atom.get_chemical_symbols()
positions=atom.get_positions()
mol_string=""
for chemical_symbol, position in zip(chemical_symbols, positions):
    mol_string+="{:s} {:.10f} {:.10f} {:.10f} \n".format(chemical_symbol, position[0], position[1], position[2])

# build a molecule
mol = gto.Mole()
mol.atom = mol_string
mol.verbose = 5
mol.output = pyscf_output
mol.unit = 'A' # angstrom
mol.charge = charge
mol.spin = spin
mol.symmetry = False

# basis set
mol.basis = basis

# define ecp
mol.ecp = ecp

# molecular build
mol.build(cart=False)  # cart = False => use spherical basis!!

# calc type setting
print(f"scf_method = {scf_method}")  # HF/DFT

if scf_method == "HF":
    # HF calculation
    if mol.spin == 0:
        print("HF kernel = RHF")
        mf = scf.RHF(mol)
        mf.chkfile = checkpoint_file
    else:
        print("HF kernel = ROHF")
        mf = scf.ROHF(mol)
        mf.chkfile = checkpoint_file

elif scf_method == "DFT":
    # DFT calculation
    if mol.spin == 0:
        print("DFT kernel = RKS")
        mf = scf.KS(mol).density_fit()
        mf.chkfile = checkpoint_file
    else:
        print("DFT kernel = ROKS")
        mf = scf.ROKS(mol)
        mf.chkfile = checkpoint_file
    mf.xc = dft_xc
else:
    raise NotImplementedError

total_energy = mf.kernel()

# HF/DFT energy
print(f"Total HF/DFT energy = {total_energy}")
print("HF/DFT calculation is done.")
print("PySCF calculation is done.")
print(f"checkpoint file = {checkpoint_file}")

You can convert the generated PySCF checkpoint file to a TREXIO file

# pyscf chkfile to TREXIO
trexio convert-from -t pyscf -i benzene.chk -b hdf5 benzene.hdf5

02 From TREXIO file to TurboRVB WF 

Next, the TREXIO file is converted to a TurboRVB wavefunction file as follows:

cd ../02_trexio_to_turborvbwf/
cp ../01_pyscf_calculation/benzene.hdf5 .

trexio-to-turborvb benzene.hdf5 -jasbasis cc-pVDZ -jascutbasis

03 JDFT ansatz - Jastrow optimization 

The next step is to optimize the Jastrow factor at the VMC level using vmcopt module. One should refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.

Copy the prepared wavefunction file fort.10 and the pseudopotential file to the work directory:

cd ../03_optimization/
cp ../02_trexio_to_turborvbwf/fort.10 fort.10
cp ../02_trexio_to_turborvbwf/pseudo.dat .
cp fort.10 fort.10_pyscf

Generate an input file for VMC optimization:

turbogenius vmcopt -g -opt_onebody -opt_twobody -opt_jas_mat -optimizer lr -vmcoptsteps 50 -steps 400 -nw 128 -reg -0.005 -num_opt_param 10

Run the VMC optimization, e,g, as follows:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius vmcopt -r

Perform the postprocess by typing:

turbogenius vmcopt -post -optwarmup 30 -plot

Check plot_energy_and_devmax.png and the files in the parameters_graphs directory.

04 JDFT ansatz - VMC 

The next step is to run a single-shot VMC calculation. This is done using the vmc module of TurboGenius. First, prepare the wavefunction and pseudopotential files:

cd ../04_vmc/
cp ../03_optimization/fort.10 fort.10
cp ../03_optimization/pseudo.dat .

Next, generate an input file datasvmc.input using:

turbogenius vmc -g -steps 3000 -nw 128

Then, run the VMC calculation, e.g., by typing:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius vmc -r

Finally, run the postprocess:

turbogenius vmc -post -bin 20 -warmup 10

Check the reblocked total energy and error in the file pip0.d.

05 JDFT ansatz - LRDMC 

Now we proceed to the lattice regularized diffusion Monte Carlo calculation that can improve a trial wavefunction obtained by a DFT calculation or a VMC optimization. One should refer to the Hydrogen tutorial for the details.

In this section, we will perform the calculation at the lattice constant alat=0.20. First, copy the prepared wavefunction and the pseudopotential files:

# LRDMC run
cd ../05_lrdmc
cp ../04_vmc/fort.10 .
cp ../04_vmc/pseudo.dat .

Next, generate an input file datasfn.input for the LRDMC calculation:

turbogenius lrdmc -g -etry -36.00 -alat -0.30 -steps 3000 -nw 128

Then, run the calculation by typing:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius lrdmc -r

Finally, run the postprocess:

turbogenius lrdmc -post -bin 20 -corr 10 -warmup 10

We will get E at a = 0.30 bohr in pip0_fn.d.

One can follow the above procedure for several choices of alat, and extrapolates the energy value at \(a \to 0\). See the Hydrogen tutorial for the concrete steps.

06 Wavefunction conversion: From JSD to JAGP 

Now we convert the wavefunction format (i.e., the degree of freedoms) from JSD to JAGP ansatz.

Copy the prepared wavefunction file fort.10 and the pseudopotential file to the work directory:

cd ../06_wavefunction_conversion_jsd_jagp/
cp ../04_vmc/fort.10 fort.10
cp ../04_vmc/pseudo.dat .

turbogenius convertwf -to agps

Warning

the original fort.10 is renamed to fort.10_bak

Please check the overlap square in out_conv:

# grep Overlap out_conv
....
Overlap square with no zero  0.9999....

Overlap square should be close to unity, i.e., if the conversion is perfect, this becomes unity.

The obtained JAGP wavefunction is fort.10.

07 JAGP ansatz - Jastrow and determinant optimization 

The next step is to optimize the Jastrow factor and the determinant level at the VMC level using vmcopt module. Here, only needed commands are shown.

Copy the prepared wavefunction file fort.10 and the pseudopotential file to the work directory:

cd ../07_optimization
cp ../06_wavefunction_conversion_jsd_jagp/fort.10 fort.10
cp ../06_wavefunction_conversion_jsd_jagp/pseudo.dat .

Generate an input file for VMC optimization:

turbogenius vmcopt -g -opt_onebody -opt_twobody -opt_jas_mat -opt_det_mat -optimizer lr -vmcoptsteps 50 -steps 400 -nw 128 -reg -0.005 -num_opt_param 10

Run the VMC optimization, e,g, as follows:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius vmcopt -r

Perform the postprocess by typing:

turbogenius vmcopt -post -optwarmup 30 -plot

Check plot_energy_and_devmax.png and the files in the parameters_graphs directory.

08 JAGP ansatz - VMC 

The next step is to run a single-shot VMC calculation. This is done using the vmc module of TurboGenius. First, prepare the wavefunction and pseudopotential files:

cd ../08_vmc/
cp ../07_optimization/fort.10 fort.10
cp ../07_optimization/pseudo.dat .

Next, generate an input file datasvmc.input using:

turbogenius vmc -g -steps 3000 -nw 128

Then, run the VMC calculation, e.g., by typing:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius vmc -r

Finally, run the postprocess:

turbogenius vmc -post -bin 20 -warmup 10

Check the reblocked total energy and error in the file pip0.d.

09 JAGP ansatz - LRDMC 

Now we proceed to the lattice regularized diffusion Monte Carlo calculation.

In this section, we will perform the calculation at the lattice constant alat=0.20. First, copy the prepared wavefunction and the pseudopotential files:

# LRDMC run
cd ../09_lrdmc
cp ../08_vmc/fort.10 .
cp ../08_vmc/pseudo.dat .

Next, generate an input file datasfn.input for the LRDMC calculation:

turbogenius lrdmc -g -etry -36.00 -alat -0.30 -steps 3000 -nw 128

Then, run the calculation by typing:

TURBOVMC_RUN_COMMAND="mpirun -np 16 turborvb-mpi.x"
export TURBOVMC_RUN_COMMAND

turbogenius lrdmc -r

Finally, run the postprocess:

turbogenius lrdmc -post -bin 20 -corr 10 -warmup 10

We will get E at a = 0.30 bohr in pip0_fn.d.

10 Summary 

The total energy obtained at the steps above and the reference value are summarized as follows:

DFT (LDA) = -36.5255 Ha
VMC (JSD) = -37.3896(28) Ha
LRDMC (JSD) = -37.6124(36) Ha (a = 0.30 bohr)
VMC (JAGP) = -37.4451(24) Ha
LRDMC (JAGP) = -37.6205(34) Ha (a = 0.30 bohr)