Quantum Package and Champ

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

First, 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_dft

These 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:

/project/project_465000321/tutorial-champ/example0*_H2O*

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: \(\exp(J)=1\).

/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.out

You 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 (\(\mathbf{r}\)) and nuclear (\(\mathbf{R}\)) coordinates. Its defined as \(\exp(J(\mathbf{r},\mathbf{R}))\), where

\[ J = f_{en} + f_{ee} + f_{een} \]

Electron-nucleus and electron-electron: \(R={1-e^{-\kappa r} \over \kappa}\)

\[ f_{en} = \sum_{i=1}^{N_{\rm elec}} \sum_{\alpha=1}^{N_{\rm nuc}} \left( {a_1 R_{i\alpha} \over 1+a_2R_{i\alpha}} + \sum_{p=2}^{N^a_{\rm ord}} a_{p+1} R_{i\alpha}^p \right) \]

\[ f_{ee} = \sum_{i=2}^{N_{\rm elec}} \sum_{j=1}^{i-1} \left( {b_1 R_{ij} \over 1+b_2R_{ij}} + \sum_{p=2}^{N^b_{\rm ord}} b_{p+1} R_{ij}^p \right) \]

Electron-electron-nucleus: \(R=\exp\left(-\kappa r \right)\)

\[ f_{een} = \sum_{i=2}^{N_{\rm elec}} \sum_{j=1}^{i-1} \sum_{\alpha=1}^{N_{\rm nuc}} \sum_{p=2}^{N^c_{\rm ord}} \sum_{k=p-1}^0 \sum_{l=l_{\rm max}}^0 c_n R_{ij}^k (R_{i\alpha}^l+R_{j\alpha}^l) (R_{i\alpha}R_{j\alpha})^m \]

where \(m={p-k-l \over 2}\)

• Typically \(N^a_{\rm ord}=N^b_{\rm ord}=5\). If \(f_{een}\) is included, \(N^c_{\rm ord}=5\).
• Dependence among \(\{c_n\}\) \(\rightarrow\) \(f_{een}\) does not contribute to cusp-conditions
• \(f_{en}\) and \(f_{een}\): different \(\{a_n\}\) and \(\{c_n\}\) for different atom types

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.

cp jastrow_optimal.1.iter20 jastrow.hf_optimal_2body

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 \(20 \times 200 / 100\) steps. Since the correlation time is less than 2 step in VMC, your walkers will be decorrelated.

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*out

In 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