@@html: <img src=”./images/QP.png” width=”100px” /> @@
We will first use Quantum Package (QP) to generate two single-determinant wave functions for the water molecule. A first one with Hartree-Fock orbitals, and a second one with PBE Kohn-Sham orbitals. Then, we will export these wave functions into the TREXIO format, which is a general format for storing arbitrary wave functions.
In a second step, we will use CHAMP to run a VMC calculation with both wave functions. We will then optimize a Jastrow factor and run DMC calculations.
For QMC calculations, we need to use pseudopotentials optimized specifically for QMC, and basis sets optimized to be used with these pseudopotentials. Here, we use the Burkatzki-Filippi-Dolg (BFD) ones except for hydrogen (the hydrogen pseudo on the website is too soft and not sufficiently accurate).
QP can read basis sets and pseudopotentials from files in GAMESS format, if the files exist in the current directory. Otherwise, it will try to look into its own database of basis sets and pseudopotentials.
Create a file called h2o.xyz: with the geometry of the water molecule:
3
Water
O 0. 0. 0.
H -0.756950272703377558 0. -0.585882234512562827
H 0.756950272703377558 0. -0.585882234512562827
Store the pseudopotential parameters in a file named PSEUDO:
H GEN 0 0
3
1.000000000000 1 25.000000000000
25.000000000000 3 10.821821902641
-8.228005709676 2 9.368618758833
O 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
Store the basis set parameters in a file named BASIS:
HYDROGEN
s 3
1 6.46417546 0.063649375945
2 1.13891461 0.339233210576
3 0.28003249 0.702654522063
s 1
1 0.05908405 1.00000000
p 1
1 0.51368060 1.00000000
OXYGEN
s 9
1 0.125346 0.055741
2 0.268022 0.304848
3 0.573098 0.453752
4 1.225429 0.295926
5 2.620277 0.019567
6 5.602818 -0.128627
7 11.980245 0.012024
8 25.616801 0.000407
9 54.775216 -0.000076
s 1
1 0.258551 1.000000
p 9
1 0.083598 0.044958
2 0.167017 0.150175
3 0.333673 0.255999
4 0.666627 0.281879
5 1.331816 0.242835
6 2.660761 0.161134
7 5.315785 0.082308
8 10.620108 0.039899
9 21.217318 0.004679
p 1
1 0.267865 1.000000
d 1
1 1.232753 1.000000
Create the EZFIO directory with the geometry, basis and pseudopotential parameters:
qp create_ezfio --pseudo=PSEUDO --basis=BASIS h2o.xyz --output=h2o_hfRun the Hartree-Fock calculation
qp_srun scf h2o_hf | tee h2o_hf.outExport the wave function into TREXIO format
qp set trexio trexio_file h2o_hf.trexio
qp_srun export_trexio h2o_hfCreate the EZFIO directory with the geometry, basis and pseudopotential parameters:
qp create_ezfio --pseudo=PSEUDO --basis=BASIS h2o.xyz --output=h2o_dftSpecify that you want to use the PBE functional.
qp set dft_keywords exchange_functional pbe
qp set dft_keywords correlation_functional pbeThe default DFT grid is very fine. We can specify we want a coarser grid to accelerate the calculations:
qp set becke_numerical_grid grid_type_sgn 1Run the Kohn-Sham calculation
qp_srun ks_scf h2o_dft | tee h2o_dft.outExport the wave function into TREXIO format
qp set trexio trexio_file h2o_dft.trexio
qp_srun export_trexio h2o_dftFirst, we can compute with QP the energies of the single-determinant wave functions with the 2 different sets of MOs.
qp_srun print_energy h2o_hf
qp_srun print_energy h2o_dftThese commands return the energy of the wavefunction contained in the EZFIO database. These values will be useful for checking that the QMC setup is OK. You should obtain the energies:
| HF MOs | -16.950384 |
| DFT MOs | -16.946588 |
Make sure that you are able to run the QP-QMC sequence, so use your own trexio file!
We advise to setup different directories for the different examples.
We will now convert the TREXIO files into input files suitable for CHAMP.
Create a new directory named H2O_HF and move the TREXIO file
h2o_hf.trexio into it. Go inside this directory and extract the
wave function information from the TREXIO file:
mkdir H2O_HF
mv h2o_hf.trexio H2O_HF
cd H2O_HF
python3 /project/project_465000321/champ/tools/trex_tools/trex2champ.py \
--trex "h2o_hf.trexio" \
--motype "Canonical" \
--backend "HDF5" \
--basis_prefix "BFD-cc-pVDZ" \
--lcao \
--geom \
--basis \
--ecp \
--det
Many files were created. Now, create a directory named pool, and
move some files into the pool:
mkdir pool
mv *.xyz *bfinfo BFD-* ECP* pool
You can now create an input file for CHAMP vmc_h2o_hf.inp :
%module general
title 'H2O HF calculation'
pool './pool/'
pseudopot ECP
basis BFD-cc-pVDZ
mode 'vmc_one_mpi1'
%endmodule
load molecule $pool/champ_v2_h2o_hf_geom.xyz
load basis_num_info $pool/champ_v3_h2o_hf_basis_pointers.bfinfo
load orbitals champ_v3_h2o_hf_trexio_orbitals.lcao
load determinants champ_v2_h2o_hf_determinants.det
load jastrow jastrow.start
%module electrons
nup 4
nelec 8
%endmodule
%module blocking_vmc
vmc_nstep 20
vmc_nblk 10000
vmc_nblkeq 1
vmc_nconf_new 0
%endmodule
Create the file for the Jastrow factor as follows, and save it as jastrow.start:
jastrow_parameter 1 0 1 0 norda,nordb,nordc 0.60000000 0.00000000 scalek,a21 0.00000000 0.00000000 (a(iparmj),iparmj=1,nparma) 0.00000000 0.00000000 (a(iparmj),iparmj=1,nparma) 0.00000000 1.00000000 (b(iparmj),iparmj=1,nparmb) (c(iparmj),iparmj=1,nparmc) (c(iparmj),iparmj=1,nparmc) end
This files implies that there is no Jastrow factor:
/project/project_465000321/tutorial-champ/example01_H2O_HF/lumi_example01.sh
is the submission script presented in section 1. Copy it in the
current dircetory and submit the job using sbatch
lumi_example01.sh. The output will be stored in vmc_h2o_hf.out
and you can grep the total energy as
grep 'total E' vmc_h2o_hf.outYou should obtain the Hartree-Fock energy.
The energies obtained with VMC without the Jastrow factor should be the same as those computed by QP at the beginning of this section.
The Jastrow factor depends on the electronic (
\[ J = fen + fee + feen \]
Electron-nucleus and electron-electron: $R={1-e-κ r \over κ}$
\[ fen = ∑i=1N_{\rm elec} ∑α=1N_{\rm nuc} \left( {a_1 Riα \over 1+a_2Riα} + ∑p=2N^a_{\rm ord} ap+1 Riα^p \right) \]
\[ fee = ∑i=2N_{\rm elec} ∑j=1i-1 \left( {b_1 Rij \over 1+b_2Rij} + ∑p=2N^b_{\rm ord} bp+1 Rij^p \right) \]
Electron-electron-nucleus:
\[ feen = ∑i=2N_{\rm elec} ∑j=1i-1 ∑α=1N_{\rm nuc} ∑p=2N^c_{\rm ord} ∑k=p-1^0 ∑l=l_{\rm max}^0 c_n Rij^k (Riα^l+Rjα^l) (RiαRjα)^m \]
where
- Typically $N^a\rm ord=N^b\rm ord=5$. If $feen$ is included, $N^c\rm ord=5$.
- Dependence among
$\{c_n\}$ $→$ $feen$ does not contribute to cusp-conditions - $fen$ and $feen
$: different $ \{a_n\}$ and$\{c_n\}$ for different atom types
- $N^a\rm ord=5$
Since we are using pseudopotentials (no e-n cusps), we always leave
$a_1=a_2=0$ and add $a_3 (riα^2), \ldots, a_6 (riα^5)$ equal to zero, which we then optimize. We do so for each atom type. - $N^b\rm ord=5$
We set
$b_1=0.5$ (for up-down e-e cusp condition), and add$b_3$ ($rij^2$),$\ldots$ ,$b_6$ ($rij^5$) equal to zero, which we then optimize.$b_1$ is modified to 0.25 for up-up and down-down electrons.The following file is your starting Jastrow factor
jastrow.start:jastrow_parameter 1 5 5 0 norda,nordb,nordc 0.60000000 scalek 0.00000000 0.00000000 0. 0. 0. 0. (a(iparmj),iparmj=1,nparma) ! e-n O 0.00000000 0.00000000 0. 0. 0. 0. (a(iparmj),iparmj=1,nparma) ! e-n H 0.50000000 1. 0. 0. 0. 0. (b(iparmj),iparmj=1,nparmb) ! e-e (c(iparmj),iparmj=1,nparmc) ! e-e-n O (c(iparmj),iparmj=1,nparmc) ! e-e-n H end
Copy your HF directory to a new directory H2O_HF_optjas2body where you optimize the Jastrow factor. Create the file jastrow.der:
jasderiv 4 4 5 0 0 0 0 nparma,nparmb,nparmc,nparmf 3 4 5 6 (iwjasa(iparm),iparm=1,nparma) ! e-n O 3 4 5 6 (iwjasa(iparm),iparm=1,nparma) ! e-n H 2 3 4 5 6 (iwjasb(iparm),iparm=1,nparmb) ! e-e 3 5 7 8 9 11 13 14 15 16 17 18 20 21 23 (c(iparmj),iparmj=1,nparmc) 3 5 7 8 9 11 13 14 15 16 17 18 20 21 23 (c(iparmj),iparmj=1,nparmc) end
where you are telling CHAMP to optimize
Now, specify the name of the info of the derivatives of the Jastrow
in the input file, below the line where the jastrow.start file is
specified. You also need to add a block with different options for the
optimizer as follows.
load jastrow jastrow.start
load jastrow_der jastrow.der
%module optwf
ioptwf 1
ioptci 0
ioptjas 1
ioptorb 0
method 'sr_n'
nopt_iter 20
nblk_max 1000
ncore 0
nextorb 100
sr_tau 0.05
sr_eps 0.001
sr_adiag 0.01
%endmodule
Optimization of the Jastrow doesn’t require a long QMC simulation in the first SR steps.
You can reduce the number of blocks in blocking_vmc to 10, and the code
will slowly increase the number of blocks to nblk_max in the optwf module.
%module blocking_vmc
vmc_nstep 20
vmc_nblk 10
vmc_nblkeq 1
vmc_nconf_new 0
%endmodule
If you grep 'total E' in the output file, you will see the optimization progressing and
generating new Jastrow factors in jastrow_optimal.1.iterX.
If you grep nblk, you will see that the code automatically increases the
maximum number of blocks.
In this section, you can use the SLURM script provided in section 1 adapted to DMC simulations. Don’t forget to edit the submission script to update the names of the input files.
Let us start to run a DMC simulation with the HF orbitals and the optimal Jastrow factor you have just generated.
Create a new directory H2O_HF_dmc2body_tau0.05 with the content of the previous HF directory with optimized Jastrow. You should now use the optimal Jastrow factor.
First, generate an input file as before where you read the wave function files (careful to load the new Jastrow factor) and perform a short VMC calculation to generate the walkers for DMC.
To shorten the VMC run, you can choose a small vmc_nblk in the main input
file and modify vmc_nconf_new to be the number of walkers per core you wish.
Here, we use the same values as for the starting iterations of the Jastrow factor
optimization:
%module blocking_vmc
vmc_nstep 20
vmc_nblk 200
vmc_nblkeq 1
vmc_nconf_new 100
%endmodule
This will generate 100 walkers per core (vmc_nconf_new) by writing
the coordinates of a walker every
A bunch of mc_configs_newX files will appear in your directory, each
containing 100 walkers.
cat mc_configs_new* >> mc_configs
rm mc_configs_new*mc_configs contains now all walkers.
Generate a DMC input
%module blocking_dmc
dmc_nstep 60
dmc_nblk 40
dmc_nblkeq 1
dmc_nconf 100
%endmodule
%module dmc
tau 0.05d0
etrial -17.24d0
icasula -1
%endmodule
You also need to change the mode keyword in the input file:
mode 'dmc_one_mpi1'
within the general module.
In the installed version of CHAMP, some debug files are being created. You can just erase them.
rm problem* walkalize*
rm mc_configs_new*To look at the energy, you can do
grep '( 100) =' dmc*outIn the last column, you have the correlation time.
Also perform another calculation with a smaller time step (e.g. 0.02).
Repeat the DMC calculation with the DFT orbitals. Compare the VMC and DMC energies.
Finally, starting from the DFT orbitals and the optimal two-body Jastrow optimize the full wave function (Jastrow and orbitals).
To this aim, set ioptorb to 1 in the optwf module.
ioptorb 1
Multiple files with geometries, basis sets and pseudopotentials can be downloaded here: Examples