02_01Lithium_dimer
00 Introduction
From this tutorial, you can learn how to calculate all-electron Variational Monte Carlo (VMC) and lattice regularized diffusion Monte Carlo (LRDMC) energies of the Li2 dimer with various ansatz. JDFT, JSD, JsAGPs, JAGPu, and JAGP(JPf).
Reference values (https://aip.scitation.org/doi/suppl/10.1063/1.3288054)
- Li atom
HF = -7.4327 Ha
Exact = -7.4781 Ha
- Li2 dimer
dbond = 2.67330 angstrom (J. Chem. Phys. 129 204105)
HF = -14.8715 Ha
Exact = -14.9951 Ha
Ebond = 24.4 kcal/mol = 38.9 mHa
01 Li2 dimer and Li atom - JDFT ansatz
01-01 Li2 dimer: Preparing a wave function
Prepare Li2 dimer structure, Li2.xyz:
2
Li2-dimer xyz file
Li 0.0000 0.0000 -1.33665000000000000000
Li 0.0000 0.0000 1.33665000000000000000
and prepare makefort10.input and convertfort10mol.input and
&system
posunits='bohr'
natoms=2
ntyp=1
pbcfort10=.false.
/
&electrons
twobody=-15
twobodypar=1.00
onebodypar=1.00
no_4body_jas=.false.
neldiff=0
/
&symmetries
eqatoms=.true.
rot_det=.true.
symmagp=.true.
/
ATOMIC_POSITIONS
3.0000000000000000 3.0000000000000000 0.0000000000000000 0.0000000000000000 -2.5259024244810400
3.0000000000000000 3.0000000000000000 0.0000000000000000 0.0000000000000000 2.5259024244810400
/
ATOM_3
&shells
nshelldet=11
nshelljas=5
/
1 1 16
1 14.24
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 16
1 0.07345
1 1 16
1 0.02805
1 1 36
1 1.534
1 1 36
1 0.2749
1 1 36
1 0.07362
1 1 36
1 0.02403
1 1 68
1 0.1144
# Parameters atomic Jastrow wf
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 36
1 0.2749
1 1 36
1 0.07362
&control
epsdgm=-1d-14
/
&mesh_info
ax=0.20
nx=64
ny=68
nz=128
/
&molec_info
nmol=3
nmolmax=3
nmolmin=3
/
You may obtain makefort10.input. Next, you can run makefort10.sh:
./makefort10.sh
which contains
makefort10.x < makefort10.input > out_make
mv fort.10_new fort.10
mv fort.10 fort.10_in
convertfort10mol.x < convertfort10mol.input > out_mol
mv fort.10_new fort.10
The generated wave function fort.10 is a symmetric Jastrow Antisymmetrized Geminal Power (JsAGPs).
01-02 Li2 dimer: generate a trial wave function using DFT.
The next step is to generate a trial wave function using the built-in DFT code. This is the minimal input file:
&simulation
itestr4=-4
iopt=1
double_mesh=.true
/
&pseudo
/
&vmc
/
&optimization
molopt=1
/
&readio
/
¶meters
/
&molecul
ax=0.2
ay=0.2
az=0.2
nx=64
ny=64
nz=128
/
&dft
maxit=50
epsdft=1d-5
mixing=1.0d0
typedft=1
optocc=0
nelocc=3
l0_at=1.0
scale_z=5
scale_hartree=-1.00
corr_hartree=.true.
linear=.false.
/
2 2 2
Then, you can run DFT by typing:
prep-serial.x < prep.input > out_prep
You can optimize one-body Jastrow.
b_onebody_list="0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5"
root_dir=`pwd`
for b_onebody in $b_onebody_list
do
mkdir onebody_$b_onebody
cd onebody_$b_onebody
cp $root_dir/template/fort.10 ./
cp $root_dir/template/prep.input ./
sed -i -e "s/b_onebody/$b_onebody/g" fort.10
prep-serial.x < prep.input > out_prep
cd $root_dir
done
echo "b_onebody dft energy" > result.out
for b_onebody in $b_onebody_list
do
cd onebody_$b_onebody
energy=`grep "Final self consistent energy" ./out_prep | awk '{print $7}'`
echo "${b_onebody} ${energy}" >> ../result.out
cd $root_dir
done
You may get
% cat result.out
b_onebody dft energy
0.5 -14.740030455469444
0.6 -14.746632760686529
0.7 -14.749842736580408
0.8 -14.750834382315727
0.9 -14.750935229881659 <- the lowest one
1.0 -14.750763468915963
1.1 -14.750537472297879
1.2 -14.750319868767475
1.3 -14.750124626530374
1.4 -14.749953844050875
1.5 -14.749807044323163
Here, the DFT double-grid integration scheme is employed.
l0_at The radius where the double-grid is used (Bohr)
scale_z The number of grids used for the double grid. Indeed, the grid sizes in the double-grid region are ax/scale_z, ay/scale_z, and ax/scale_z
scale_hartree, corr_hartree, and linear can be default values.
01-03 Li2 dimer: Jastrow factor optimization (WF=JDFT)
One should refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory and copy the trial wavefunction generated by the DFT calculation at the optimal onebody Jastrow term:
cd 02optimization cp ../01trial_wavefunction/02onebody_jastrow_opt/onebody_0.9/fort.10_new fort.10 cp fort.10 fort.10_dft
Prepare an input file
datasmin.inputStart a VMC optimization:
turborvb-serial.x < datasmin.input > out_min
After finishing the calculation, you can delete template files:
rm -r turborvb.scratch/
Confirm energy convergence by typing:
plot_Energy.sh out_minCheck the convergence of devmax by typing:
plot_devmax.sh out_minNext step is to average optimized variational parameters.
readalles.x << ____EOS 1 41 1 0 0 1000 ____EOS
01-04 Li2 dimer: VMC (WF=JDFT)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory, and copy the trial wavefunction:
cd 03vmc/ cp ../02optimization/fort.10 ./
Then, copy
datasmin.inputand rename it asave.in:cp ../02optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 20 5 1
A reblocked total energy and its variance is written in
pip0.d.% cat pip0.d ... number of bins read = 14996 Energy = -14.9757876300445 4.853630781726438E-004 Variance square = 7.979609547808308E-002 2.379635260709342E-004 Est. energy error bar = 4.748991835467470E-004 2.848268581908139E-006 Est. corr. time = 3.39078657629560 3.817021016290687E-002
VMC (JDFT) = -14.9759(9) Ha
01-05 Li2 dimer: LRDMC (WF=JDFT)
One should refer to the Hydrogen tutorial for the details. Here, only needed explanations and commands are shown.
First, prepare an input file datasfn.input for LRDMC calculation:
&simulation
itestr4=-6
ngen=200000
iopt=2
maxtime=345600
/
&pseudo
/
&dmclrdmc
tbra=0.1d0
etry=-15.00d0
alat=0.XXd0
!alat2=0.0d0
iesrandoma=.true.
/
&readio
/
¶meters
/
Here are brief explanations of the variables for the LRDMC calculation:
&dmclrdmc section
tbra projection time (i.e, \(\exp(-\tau \cdot \hat{\mathcal{H}})\)). Set 0.1 in general. However, for a heavy element, it is better to choose a smaller value. Please check Average number of survived walkers in out_fn
Av. num. of survived walkers/ # walkers in the branching
0.9939
If the number is too small, try smaller tbra.
etry Put a DFT of VMC energy. \(\Gamma\) in Eq.6 of the review paper is set 2 \(\times\) etry
alat The lattice space for discretizing the Hamitonian. If you do a single grid calculation (i.e., alat2=0.0d0), please put a negative value. If you do a double-grid calculation (See. Nakano’s paper), put a positive value and set iesrandoma=.true.. This trick is needed for satisfying the detailed-valance condition.
alat2 The corser lattice space used in the double-grid calculation. If you put 0.0d0, TurboRVB does a single grid calculation. If you want to do a double-grid calculation for a compound include Z > 2 element, please comment out alat2 because alat2 is automatically set (See Nakano’s paper).
Prepare different working directories, copy fort.10 to each directory, and set the corresponding alat.
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
And then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd 03lrdmc/
cp ../02vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
One has collected all LRDMC energis into evsa.in
# Z=3 (Li)
2 4 4 1
0.22222 -14.9907899448815 3.101461151525789E-004 <- 1/(1.5*Z)
0.27778 -14.9920193657668 3.175522770991489E-004 <- 1/(1.2*Z)
0.33333 -14.9933364338623 3.387085675247969E-004 <- 1/(1.0*Z)
0.41667 -14.9969823559963 3.650299072111256E-004 <- 1/(0.8*Z)
funvsa.x is a tool for a quadratic fitting:
funvsa.x < evsa.in > evsa.out
You can see
Reduced chi^2 = 0.258290367368527
Coefficient found
1 -14.9894036676827 1.126607950073031E-003 <- E_0
2 -2.333404434142712E-002 2.310746587271902E-002 <- k_1
3 -0.116381701497219 0.102032433967512 <- k_2
- Li2 dimer
HF = -14.8715 Ha
VMC(JDFT) = -14.9759(9) Ha
LRDMC(JDFT) = -14.9894(11) Ha
Exact = -14.9951 Ha
01-06 Li atom: preparation of a wave function
First, prepare makefort10.input and convertfort10mol.input:
makefort10.input:
&system
posunits='bohr'
natoms=1
ntyp=1
pbcfort10=.false.
/
&electrons
twobody=-15
twobodypar=1.00
onebodypar=0.90
no_4body_jas=.false.
neldiff=1 ! the number of unpaired electrons !!
/
&symmetries
eqatoms=.true.
rot_det=.true.
symmagp=.true.
/
ATOMIC_POSITIONS
3.0000000000000000 3.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
/
ATOM_3
&shells
nshelldet=11
nshelljas=5
/
1 1 16
1 14.24
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 16
1 0.07345
1 1 16
1 0.02805
1 1 36
1 1.534
1 1 36
1 0.2749
1 1 36
1 0.07362
1 1 36
1 0.02403
1 1 68
1 0.1144
# Parameters atomic Jastrow wf
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 36
1 0.2749
1 1 36
1 0.07362
convertfort10mol.input:
&control
epsdgm=-1d-14
/
&mesh_info
ax=0.20
nx=64
ny=64
nz=64
/
&molec_info
nmol=1
nmolmax=1
nmolmin=1
/
Next, you can run makefort10.sh:
% ./makefort10.sh
which contains
makefort10.x < makefort10.input > out_make
mv fort.10_new fort.10
mv fort.10 fort.10_in
convertfort10mol.x < convertfort10mol.input > out_mol
mv fort.10_new fort.10
The generated wave function fort.10 is a symmetric Jastrow Antisymmetrized Geminal Power (JsAGPs).
01-07 Li atom: generate a trial wave function using DFT.
One should refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
prep-serial.x < prep.input > out_prep
01-08 Li atom: Jastrow factor optimization (WF=JDFT)
One should refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory and copy the trial wavefunction generated by the DFT calculation at the optimal onebody Jastrow term:
cd 02optimization cp ../01trial_wavefunction/01DFT/fort.10_new fort.10 cp fort.10 fort.10_dft
Prepare an input file
datasmin.inputStart a VMC optimization:
turborvb-serial.x < datasmin.input > out_min
After finishing the calculation, you can delete template files:
rm -r turborvb.scratch/
Confirm energy convergence by typing:
plot_Energy.sh out_minCheck the convergence of devmax by typing:
plot_devmax.sh out_minNext step is to average optimized variational parameters.
readalles.x << ____EOS 1 41 1 0 0 1000 ____EOS
01-09 Li atom: VMC (WF=JDFT)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory, and copy the trial wavefunction:
cd 03vmc/ cp ../02optimization/fort.10 ./
Then, copy
datasmin.inputand rename it asave.in:cp ../02optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 20 5 1
A reblocked total energy and its variance is written in
pip0.d.number of bins read = 14996 Energy = -7.47539042703452 2.682804200575495E-004 Variance square = 2.847683107759093E-002 1.295272783812896E-004 Est. energy error bar = 2.683664075931440E-004 1.733280980510503E-006 Est. corr. time = 3.03421296313493 3.544802996972372E-002
VMC (JDFT) = -14.9759(9) Ha (Li2 dimer) VMC (JDFT) = -7.4754(3) Ha (Li atom)
- Li2 dimer
Ebond = 38.9 mHa (experiment)
Ebond = 25.1 mHa (VMC-JDFT)
01-10 Li atom: LRDMC (WF=JDFT)
One should refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Prepare different working directories for each value of alat, copy fort.10 to each directory, and create an input file datasfn.input with the corresponding alat:
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
Then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd 03lrdmc/
cp ../02vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
One has collected all LRDMC energis into evsa.in
# Z=3 (Li)
2 4 4 1
0.22222 -7.47824129633644 1.450042242028465E-004 <- 1/(1.5*Z)
0.27778 -7.47833093539214 1.562945930953174E-004 <- 1/(1.2*Z)
0.33333 -7.47909688312053 1.628394384627272E-004 <- 1/(1.0*Z)
0.41667 -7.48032309244573 1.866069140147947E-004 <- 1/(0.8*Z)
funvsa.x is a tool for a quadratic fitting:
funvsa.x < evsa.in > evsa.out
You can see
Reduced chi^2 = 1.52438473425829
Coefficient found
1 -7.47794022870940 5.344739888993537E-004 <- E_0
2 -1.467218812266460E-003 1.107590136996555E-002 <- k_1
3 -7.145382242986854E-002 4.939621583555785E-002 <- k_2
- Li2 dimer
LRDMC(JDFT) = -14.9894(11) Ha
LRDMC(JDFT) = -7.4779(5) Ha
Ebond = 38.9 mHa (experiment)
Ebond = 25.1 mHa (VMC-JDFT)
Ebond = 33.6 mHa (LRDMC-JDFT)
02 Li2 dimer and Li atom - JSD ansatz
02-01 Li2 dimer and Li atom: VMC-optimization (WF=JSD)
One can optimize the determinant part of the JDFT ansatz, the resultant
ansatz is called JSD.
In this case, one does not have to convert the ansatz, just put the
following section in datasmin.input at the optimization step.
datasmin.input for the Li dimer:
&optimization
molopt=-1
...
/
¶meters
iesd=1
iesfree=1
iessw=1
/
&molecul
ax=0.10
ay=0.10
az=0.10
nx=150
ny=150
nz=250
nmolmin=3
nmolmax=3
/
datasmin.input for the Li atom:
&optimization
molopt=-1
...
/
¶meters
iesd=1
iesfree=1
iessw=1
/
&molecul
ax=0.10
ay=0.10
az=0.10
nx=150
ny=150
nz=150
nmolmin=1
nmolmax=1
/
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Li2 dimer
Copy
fort.10from JDFT ansatz:cp ../../../01JDFT/01Li2_dimer/03vmc/fort.10 . cp fort.10 fort.10_jdft
Run optimization:
turborvb-serial.x < datasmin.input > out_min
Average optimized variational parameters:
readalles.x << ____EOS 1 41 1 0 0 1000 ____EOS
Li atom
Copy
fort.10from JDFT ansatz:cp ../../../01JDFT/02Li_atom/03vmc/fort.10 . cp fort.10 fort.10_jdft
Run optimization:
turborvb-serial.x < datasmin.input > out_min
Average optimized variational parameters:
readalles.x << ____EOS 1 41 1 0 0 1000 ____EOS
02-02 Li2 dimer and Li atom: VMC (WF=JSD)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory, and copy the trial wavefunction:
cd 02vmc/ cp ../01optimization/fort.10 .
Then, copy
datasmin.inputand rename it asave.in:cp ../01optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 10 2 1
A reblocked total energy and its variance is written in
pip0.d.
02-03 Li2 dimer and Li atom: LRDMC (WF=JSD)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Prepare different working directories for each value of alat, copy fort.10 to each directory, and create an input file datasfn.input with the corresponding alat:
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
Then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd ../03lrdmc/
cp ../02vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
funvsa.x < evsa.in > evsa.out
Finally we got:
Li2 dimer
HF = -14.8715 Ha
VMC(JDFT) = -14.9759(9) Ha
LRDMC(JDFT) = -14.9894(11) Ha
VMC(JSD) = -14.9803(5) Ha
LRDMC(JSD) = -14.9909(12) Ha
Exact = -14.9951 Ha
Li atom
HF = -7.4327 Ha
VMC(JDFT) = -7.4754(3) Ha
LRDMC(JDFT) = -7.4779(5) Ha
VMC(JSD) = -7.4769(2) Ha
LRDMC(JSD) = -7.4779(4) Ha
Exact = -7.4781 Ha
Binding energy
Ebond = 38.9 mHa (experiment)
03 Li2 dimer and Li atom - JsAGPs ansatz
The procedure is the almost same as in the Hydrogen-dimer tutorial.
Three hybrid-orbitals (nhyb=2) were employed here.
03-01 Li2 dimer and Li atom: Conversion of WF (WF=JsAGPs)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Li2 dimer
create a work directory, and copy
fort.10from JSD ansatz.cd 03JsAGPs/01Li2_dimer/01convert_WF_JSD_to_JAGP cp ../../../02JSD/01Li2_dimer/02vmc/fort.10 ./fort.10_in
copy
makefort10.inputfrom JDFT ansatz, and edit:cp ../../../01JDFT/01Li2_dimer/01trial_wavefunction/00makefort10/makefort10.input . set `grep -A1 "Parameters Jastrow two body" ./fort.10_in | tail -1` twobodyjas=$2 onebodyjas=$3 sed -i \ -e '/twobodypar/s/=.*$/'$twobodyjas'/' \ -e '/onebodypar/s/=.*$/'$onebodyjas'/' \ -e '/nshelljas/a ndet_hyb=2' \ makefort10.input
run makefort10
makefort10.x < makefort10.input > out_make mv fort.10_new fort.10_out
convert wavefunction
convertfort10-serial.x < convertfort10.input > out_conv grep "Overlap square Geminal" out_conv
copy Jastrow factor to wavefunction
mv fort.10_new fort.10 cp fort.10_in fort.10_new copyjas.x > out_copyjas cleanfort10.x > out_cleanfort10 cp fort.10_clean fort.10
Li atom
create a work directory, and copy
fort.10from JSD ansatz.cd 03JsAGPs/02Li_atom/01convert_WF_JSD_to_JAGP cp ../../../02JSD/02Li_atom/02vmc/fort.10 ./fort.10_in
copy
makefort10.inputfrom JDFT ansatz, and edit:cp ../../../01JDFT/02Li_atom/01trial_wavefunction/00makefort10/makefort10.input . set `grep -A1 "Parameters Jastrow two body" ./fort.10_in | tail -1` twobodyjas=$2 onebodyjas=$3 sed -i \ -e '/twobodypar/s/=.*$/'$twobodyjas'/' \ -e '/onebodypar/s/=.*$/'$onebodyjas'/' \ -e '/nshelljas/a ndet_hyb=2' \ makefort10.input
run makefort10
makefort10.x < makefort10.input > out_make mv fort.10_new fort.10_out
convert wavefunction
convertfort10-serial.x < convertfort10.input > out_conv grep "Overlap square Geminal" out_conv
copy Jastrow factor to wavefunction
mv fort.10_new fort.10 cp fort.10_in fort.10_new copyjas.x > out_copyjas cleanfort10.x > out_cleanfort10 cp fort.10_clean fort.10
03-02 Li2 dimer and Li atom: conversion check (WF=JsAGPs)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
cd ../02conversion_check/
cp ../01convert_WF_JSD_to_JAGP/fort.10 .
cp ../01convert_WF_JSD_to_JAGP/fort.10_in fort.10_corr
turborvb-serial.x < datasvmc.input > out_vmc;
readforward-serial.x < datasvmc.input > out_read;
cat corrsampling.dat
03-03 Li2 dimer and Li atom: VMC-optimization (WF=JsAGPs)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Warning
If you want to optimized the contraction coefficients, you should put
iesup=1 in the ¶meters section. The code optimizes also
contraction coefficients when you put iesup=1 with itestr4=-4
or itestr4=-9. When you put iesup=1 with itestr4=-8
or itestr4=-5, the code optimizes not only contraction coefficients
but also exponents of basis set.
create a work directory, and copy
fort.10:cd ../03optimization cp ../01convert_WF_JSD_to_JAGP/fort.10 .
prepare
datasmin.input... ¶meters iesd=1 iesfree=1 iessw=1 iesup=1 ...
run optimization
turborvb-serial.x < datasmin.input > out_min readalles.x << ____EOS 1 81 1 0 0 1000 ____EOS
03-04 Li2 dimer and Li atom: VMC (WF=JsAGPs)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory, and copy the trial wavefunction:
cd ../04vmc/ cp ../03optimization/fort.10 .
Then, copy
datasmin.inputand rename it asave.in:cp ../03optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 10 2 1
A reblocked total energy and its variance is written in
pip0.d.
03-05 Li2 dimer and Li atom: LRDMC (WF=JsAGPs)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Prepare different working directories for each value of alat, copy fort.10 to each directory, and create an input file datasfn.input with the corresponding alat:
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
Then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd ../05lrdmc/
cp ../04vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
funvsa.x < evsa.in > evsa.out
Finally we got:
Li2 dimer
HF = -14.8715 Ha
VMC(JDFT) = -14.9759(9) Ha
LRDMC(JDFT) = -14.9894(11) Ha
VMC(JSD) = -14.9803(5) Ha
LRDMC(JSD) = -14.9909(12) Ha
VMC(JsAGPs) = -14.9821(4) Ha
LRDMC(JsAGPs) = -14.991(1) Ha
Exact = -14.9951 Ha
Li atom
HF = -7.4327 Ha
VMC(JDFT) = -7.4754(3) Ha
LRDMC(JDFT) = -7.4779(5) Ha
VMC(JSD) = -7.4769(2) Ha
LRDMC(JSD) = -7.4779(4) Ha
VMC(JsAGPs) = -7.477(2) Ha
LRDMC(JsAGPs) = -7.4783(4) Ha
Exact = -7.4781 Ha
Binding energy
Ebond = 38.9 mHa (experiment)
Warning
For a real run (i.e., for a peer-reviewed paper), one should optimize variational parameters much more carefully. We recommend that one consult to an expert or a developer of TurboRVB, or carefully read the Wavefuntion optimization part.
04 Li2 dimer and Li atom - JAGPu ansatz
04-01 Li2 dimer: prepration of a wave function
As in the JDFT procedure, prepare Li2 dimer structure,
Li2.xyz :
2
Li2-dimer xyz file
Li 0.0000 0.0000 -1.33665000000000000000
Li 0.0000 0.0000 1.33665000000000000000
The point is that you should set twobody=-22 and symmagp=.false. in makefort10.input:
&system
posunits='bohr'
natoms=2
ntyp=1
pbcfort10=.false.
/
&electrons
twobody=-22
twobodypar=1.00
onebodypar=0.90
no_4body_jas=.false.
neldiff=0
/
&symmetries
nosym=.true.
eqatoms=.true.
!rot_det=.true.
symmagp=.false.
/
ATOMIC_POSITIONS
3.0000000000000000 3.0000000000000000 0.0000000000000000 0.0000000000000000 -2.5259024244810400
3.0000000000000000 3.0000000000000000 0.0000000000000000 0.0000000000000000 2.5259024244810400
/
ATOM_3
&shells
nshelldet=11
nshelljas=5
/
1 1 16
1 14.24
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 16
1 0.07345
1 1 16
1 0.02805
1 1 36
1 1.534
1 1 36
1 0.2749
1 1 36
1 0.07362
1 1 36
1 0.02403
1 1 68
1 0.1144
# Parameters atomic Jastrow wf
1 1 16
1 4.581
1 1 16
1 1.58
1 1 16
1 0.564
1 1 36
1 0.2749
1 1 36
1 0.07362
convertfort10mol.input reads as follows:
&control
epsdgm=-1d-14
/
&mesh_info
ax=0.20
nx=64
ny=64
nz=128
/
&molec_info
nmol=3
nmolmax=3
nmolmin=3
/
The other procedure is the same. You may obtain makefort10.input. Next, you can run makefort10.sh:
% ./makefort10.sh
which contains
makefort10.x < makefort10.input > out_make
mv fort.10_new fort.10
mv fort.10 fort.10_in
convertfort10mol-serial.x < convertfort10mol.input > out_mol
mv fort.10_new fort.10
The generated wave function fort.10 is a Jastrow Antisymmetrized Geminal Power including a triplet correlation (JAGPu).
Note
If you want to use a more general spin-dependent Jastrow, i.e., independent parameters for the parallel and opposite spins (twobody = -27), you should manually put variational parameters for the two-body Jastrow part after generating fort.10. Since makefort10.x supports only one variable for the two-body part at present, while -27 has two independent two-body Jastrow variational parameters.
Indeed, If you want to use a more general spin-dependent Jastrow, whenever you generate a new fort.10 file, you should replace
# Parameters Jastrow two body
-2 1.00000000000000 0.900000000000000
in a generated fort.10 with
# Parameters Jastrow two body
-3 1.00000000000000 1.00000000000000 0.900000000000000
If the number of species in a system is more than 1, e.g., LiH, you should also put independent one-body Jastrows. Namely, you should replace
# Parameters Jastrow two body
-2 1.00000000000000 0.900000000000000
in a generated fort.10 with
# Parameters Jastrow two body
-4 1.00000000000000 1.00000000000000 0.900000000000000 0.900000000000000
04-02 Li2 dimer: generate a trial wave function using DFT.
The next step is to generate a trial wave function using the built-in DFT code. This is the minimal input file:
&simulation
itestr4=-4
iopt=1
double_mesh=.true.
/
&pseudo
/
&vmc
/
&optimization
molopt=1
/
&readio
/
¶meters
/
&molecul
ax=0.2
ay=0.2
az=0.2
nx=64
ny=64
nz=128
/
&dft
maxit=50
epsdft=1d-5
mixing=1.0d0
typedft=4
optocc=1
l0_at=1.0
scale_z=5
scale_hartree=-1.00
corr_hartree=.true.
linear=.false.
nxs=1
nys=1
nzs=2
randspin=-1.0d0
h_field=1.0d-2
/
1 -1
Then, you can run DFT by typing:
prep-serial.x < prep.input > out_prep
04-03 Li2 dimer: Check local magnetic moments
You can check the obtained magnetic moments by using plot_orbitals.x.
First, copy fort.10 obtained by DFT calculation:
cp ../01DFT/fort.10_new fort.10
Then, plot spin density:
% plot_orbitals.x
Number of molecular orbitals : 6
Choose box size (x,y,z)
15 15 15
15.0000000000000 15.0000000000000 15.0000000000000
Choose number of mesh points (x,y,z) :
61 61 61
61 61 61
Choose orbitals to tabulate (possible answers all, partial, charge, spin) :
spin
spin
Please give the lowest molecular orbital within 1 and 6
1
Number of fully occupied molecular orbital/total number occupied by up and down?
3 3
Momentum magnetization ? (unit 2pi/cellscale)
0 0 0
You may obtain –xsf output_spin000000.xsf. Depict it using xcrysden:
% xcrysden --xsf output_spin000000.xsf
Thus, one can obtain an AFM trial wavefunction.
04-04 Li2 dimer and Li atom: Convert JDFT WF to JAGPu one
Next step is to convert the optimized JDFT ansatz to a JAGPu one. You should check the consistency after conversion.
copy
fort.10obtained from DFT calculation:cp ../01trial_wavefunction/01DFT/fort.10_new ./fort.10_in
copy
makefort10.input, and edit:cp ../01trial_wavefunction/00makefort10/makefort10.input . sed -i -e '/&symmetries/a nosym_contr=.true.' makefort10.input sed -i -e '/nshelljas/a ndet_hyb=2' makefort10.input
run makefort10
makefort10.x < makefort10.input > out_make mv fort.10_new fort.10_out
convert wavefunction
convertfort10-serial.x < convertfort10.input > out_conv mv fort.10_new fort.10
run correlated sampling for conversion check
cd ../03conversion_check cp ../02convert_WF_JDFT_to_JAGP/fort.10 . cp ../02convert_WF_JDFT_to_JAGP/fort.10_in fort.10_corr turborvb-serial.x < datasvmc.input > out_vmc readforward-serial.x < datasvmc.input > out_read
Li2 dimer:
% cat corrsampling.dat
Number of bins 146
reference energy: E(fort.10) -0.149056547E+02 0.903277814E-02
reweighted energy: E(fort.10_corr) -0.149153980E+02 0.890288018E-02
reweighted difference: E(fort.10)-E(fort.10_corr) 0.974326537E-02 0.795833755E-03
Overlap square : (fort.10,fort.10_corr) 0.993230204E+00 0.284523743E-03
Li atom:
% cat corrsampling.dat
Number of bins 146
reference energy: E(fort.10) -0.148896680E+02 0.929468673E-02
reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468634E-02
reweighted difference: E(fort.10)-E(fort.10_corr) 0.124625377E-07 0.316227766E-07
Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.316227766E-07
04-05 Li2 dimer and Li atom: Optimization (WF=JAGPu)
Now that you have obtained a good trial JAGPu wavefunction, you can optimize its nodal surface at the VMC level.
Then, copy the converted WF.
cp ../02convert_WF_JDFT_to_JAGP/fort.10 ./
cp fort.10 fort.10_dft
Prepare datasmin.input:
&simulation
itestr4=-4
ngen=200000
iopt=1
maxtime=10800
/
&pseudo
/
&vmc
/
&optimization
nweight=2000
nbinr=20
iboot=0
tpar=3.5d-1
parr=5.0d-3
npbra=5
parcutpar=4.5d0
/
&readio
/
¶meters
iesd=1
iesfree=1
iessw=1
iesup=1
/
The difference from datasmin.input in 03JsAGPs/03optimization is iessw,
because we optimize the determinant part (nodal surface) at this step.
Run a VMC-opt run.
#run vmc-opt
turborvb-serial.x < datasmin.input > out_min
#average fort.10
readalles.x << ____EOS
1 81 1 0
0
1000
____EOS
The rest of the procesure (e.g., averaging the variational parameters) is the same.
Warning
For a real run (i.e., for a peer-reviewed paper), one should optimize variational parameters much more carefully. We recommend that one consult to an expert user or a developer of TurboRVB, or carefully read the Wavefuntion optimization part.
04-06 Li2 dimer and Li atom: VMC and LRDMC
VMC and LRDMC procesures are the same as in the JsAGPs case.
VMC
First, create a working directory, and copy the trial wavefunction:
cd ../05vmc/ cp ../04ptimization/fort.10 .
Then, copy
datasmin.inputand rename it asave.in:cp ../04optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 10 2 1
A reblocked total energy and its variance is written in
pip0.d.
LRDMC
Prepare different working directories for each value of alat, copy fort.10 to each directory, and create an input file datasfn.input with the corresponding alat:
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
Then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd ../06lrdmc/
cp ../05vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
funvsa.x < evsa.in > evsa.out
Li2 dimer
HF = -14.8715 Ha
VMC(JDFT) = -14.9759(9) Ha
LRDMC(JDFT) = -14.9894(11) Ha
VMC(JsAGPs) = -14.9821(4) Ha
LRDMC(JsAGPs) = -14.991(1) Ha
VMC(JAGPu) = -14.9819(5) Ha
LRDMC(JAGPu) = -14.990(1) Ha
Exact = -14.9951 Ha
Li atom
HF = -7.4327 Ha
VMC(JDFT) = -7.4754(3) Ha
LRDMC(JDFT) = -7.4779(5) Ha
VMC(JsAGPs) = -7.477(2) Ha
LRDMC(JsAGPs) = -7.4783(4) Ha
VMC(JAGPu) = -7.477(3) Ha
LRDMC(JAGPu) = -7.4791(4) Ha
Exact = -7.4781 Ha
Binding energy
Ebond = 38.9 mHa (experiment)
05 Li2 dimer and Li atom - JAGP (JPf) ansatz
05-01 (spin-unpolarized case) Li2 dimer: Convert JDFT WF to JAGP one
The most important procedure in a Pfaffian calculation is to convert a JDFT or JAGPu ansatz to JAGP(JPf) ansatz. Since the JAGP ansatz is a special case of the JPf one, where only \(G_{ud}\) and \(G_{du}\) terms are defined as described in the section review paper, the conversion can be realized just by direct substitution. Therefore, the main challenge is to find a reasonable initialization for the two spin-triplet sectors \(G_{uu}\) and \(G_{dd}\) that are not described in the JAGP and that otherwise have to be set to 0. There are two possible approaches to convert an ansatz: (i) for polarized systems, we can build the \(G_{uu}\) block of the matrix by using an even number of \(\{ \phi_i\}\) and build an antisymmetric \(g_{uu}\), where the eigenvalues \(\lambda_k\) are chosen to be large enough to occupy certainly these unpaired states, as in the standard Slater determinant used for our initialization. Again, we emphasize that this works only for polarized systems. (ii) The second approach that also works in a spin-unpolarized case is to determine a standard broken symmetry single determinant ansatz (e.g., by TurboRVB built-in DFT within the LSDA) and modify it with a global spin rotation. Indeed, in the presence of finite local magnetic moments, it is often convenient to rotate the spin moments of the WF in a direction perpendicular to the spin quantization axis chosen for our spin-dependent Jastrow factor, i.e., the \(z\) quantization axis. In this way one can obtain reasonable initializations for \(G_{uu}\) and \(G_{dd}\). TurboRVB allows every possible rotation, including an arbitrary small one close to the identity. A particularly important case is when a rotation of \(\pi/2\) is applied around the \(y\) direction. This operation maps \(|\uparrow \rangle \rightarrow \frac{1} {\sqrt{2}} \left( | \uparrow \rangle + |\downarrow \rangle \right) \mbox{ and } |\downarrow \rangle \rightarrow \frac 1 {\sqrt{2}} \left( | \uparrow \rangle - |\downarrow \rangle \right).\) One can convert from a AGP the pairing function that is obtained from a VMC optimization \({g_{ud}}(\mathbf{i},\mathbf{j}) = {f_S}({{\mathbf{r}}_i},{{\mathbf{r}}_j})\frac{{\left| { \uparrow \downarrow } \right\rangle - \left| { \downarrow \uparrow } \right\rangle }}{{\sqrt 2 }} + {f_T}({{\mathbf{r}}_i},{{\mathbf{r}}_j})\frac{{\left| { \uparrow \downarrow } \right\rangle + \left| { \downarrow \uparrow } \right\rangle }}{{\sqrt 2 }}\) to a Pf one \({g_{ud}}(\mathbf{i},\mathbf{j}) \to g\left( {\mathbf{i},\mathbf{j}} \right){\text{ }} = {f_S}({{\mathbf{r}}_i},{{\mathbf{r}}_j})\frac{{\left| { \uparrow \downarrow } \right\rangle - \left| { \downarrow \uparrow } \right\rangle }}{{\sqrt 2 }} + {f_T}({{\mathbf{r}}_i},{{\mathbf{r}}_j})\left( {\left| { \uparrow \uparrow } \right\rangle - \left| { \downarrow \downarrow } \right\rangle } \right).\) This transformation provides a meaningful initialization to the Pfaffian WF that can be then optimized for reaching the best possible description of the ground state within this ansatz.
The strategy (ii) is employed for the Li dimer (i.e., unpolarized case) while (i) is employed for the Li atom (i.e., polarized case)
First, copy
fort.10from JAGPu ansatz, and run VMC for corrlated sampling.cp ../../../04JAGPu/01Li2_dimer/01trial_wavefunction/01DFT/fort.10_new fort.10_in cp fort.10_in fort.10_dft cp fort.10_dft fort.10 turborvb-serial.x < datasvmc.input > out_vmc rm fort.10
Convert wavefunction from JDFT ansatz to uncontracted JAGP ansatz:
cp ../../../04JAGPu/01Li2_dimer/01trial_wavefunction/00makefort10/makefort10.input makefort10_agp_uncont.input makefort10.x < makefort10_agp_uncont.input > out_make_agp_uncont mv fort.10_new fort.10_out convertfort10-serial.x < convertfort10.input > out_conv_agp_uncont mv fort.10_new fort.10_agp_uncont
Run correlated sampling
cp fort.10_dft fort.10 cp fort.10_agp_uncont fort.10_corr readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_dft_agp_uncont.dat rm fort.10
% cat corrsampling_dft_agp_uncont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468647E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.122986368E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.316227766E-07
Convert wavefunction from uncontracted JAGP ansatz to uncontracted JPf (norotate) ansatz:
cp fort.10_agp_uncont fort.10_in cp makefort10_agp_uncont.input makefort10_pf_uncont.input
Edit
makefort10_pf_uncont.inputsed -i -e '/&system/a yes_pfaff=.true.' -e '/&symetries/a rot_pfaff=.false.' makefort10_pf_uncont.input
Convert wavefunction
makefort10.x < makefort10_pf_uncont.input > out_make_pfaff_uncont mv fort.10_new fort.10_out convertpfaff.x norotate > out_conv_pfaff_norotate cp fort.10_new fort.10_pfaff_uncont
Run correlated sampling
cp fort.10_dft fort.10 cp fort.10_pfaff_uncont fort.10_corr readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_agp_pfaff_uncont.dat
% cat corrsampling_agp_pfaff_uncont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468647E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.122986368E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.316227766E-07
Convert wavefunction from uncontracted JPf ansatz to contracted JPf ansatz:
cp fort.10_pfaff_uncont fort.10_in cp makefort10_pf_uncont.input makefort10_pf_cont.input
Edit makefort10_pf_cont.input
sed -i -e '/&symmetries/a nosym_contr=.true.' -e '/nshelljas/a ndet_hyb=2' makefort10_pf_cont.input
Convert wavefunction
makefort10.x < makefort10_pf_cont.input > out_make_pfaff_cont mv fort.10_new fort.10_out convertfort10-serial.x < convertfort10.input > out_conv_pfaff_cont cp fort.10_new fort.10_pfaff_cont cp fort.10_pfaff_cont fort.10_corr
Run correlated sampling
cp fort.10_dft fort.10 readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_dft_pfaff_cont.dat
% cat corrsampling_dft_pfaff_cont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468645E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.115194734E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.505322962E-07
Rotate wavefunction by 1/8 pi:
cp fort.10_pfaff_cont fort.10_in cp fort.10_pfaff_cont fort.10_out echo "-0.125" | convertpfaff.x > out_conv_pfaff_rot cp fort.10_new fort.10_pfaff_cont_rot
Clean fort.10
cp fort.10_pfaff_cont_rot fort.10 cleanfort10.x mv fort.10_clean fort.10 cp fort.10 fort.10_final cp fort.10 fort.10_corr
Run correlated sampling
cp fort.10_dft fort.10 readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_dft_pfaff_cont_rot.dat rm fort.10
% cat corrsampling_dft_pfaff_cont_rot.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148904372E+02 0.936290264E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.769116534E-03 0.396919093E-03 Overlap square : (fort.10,fort.10_corr) 0.998728350E+00 0.512541253E-03
The total energy was a bit lost by the rotation.
Note
If you want to use a more general spin-dependent Jastrow, i.e., independent parameters for the parallel and opposite spins (twobody = -27), you should manually put variational parameters for the two-body Jastrow part whenever generating a new fort.10.
05-02 (spin-polarized case) Li atom: Convert JDFT WF to JAGP one
First, copy
fort.10from JAGPu ansatz, and run VMC for correlated sampling.cp ../../../04JAGPu/02Li_atom/01trial_wavefunction/01DFT/fort.10_new fort.10_in cp fort.10_in fort.10_dft cp fort.10_dft fort.10 turborvb-serial.x < datasvmc.input > out_vmc rm fort.10
Convert wavefunction from JDFT ansatz to uncontracted JAGP ansatz:
cp ../../../04JAGPu/02Li_atom/01trial_wavefunction/00makefort10/makefort10.input makefort10_agp_uncont.input makefort10.x < makefort10_agp_uncont.input > out_make_agp_uncont mv fort.10_new fort.10_out convertfort10-serial.x < convertfort10.input > out_conv_agp_uncont mv fort.10_new fort.10_agp_uncont
Run correlated sampling
cp fort.10_dft fort.10 cp fort.10_agp_uncont fort.10_corr readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_dft_agp_uncont.dat rm fort.10
% cat corrsampling_dft_agp_uncont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468647E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.122986368E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.316227766E-07
Convert wavefunction from uncontracted JAGP ansatz to uncontracted JPf (norotate) ansatz:
cp fort.10_agp_uncont fort.10_in cp makefort10_agp_uncont.input makefort10_pf_uncont.input
Edit makefort10_pf_uncont.input
sed -i -e '/&system/a yes_pfaff=.true.' -e '/&symetries/a rot_pfaff=.false.' makefort10_pf_uncont.input
Convert wavefunction
makefort10.x < makefort10_pf_uncont.input > out_make_pfaff_uncont mv fort.10_new fort.10_out echo "10000" | convertpfaff.x norotate > out_conv_pfaff_norotate cp fort.10_new fort.10_pfaff_uncont
Run correlated sampling
cp fort.10_dft fort.10 cp fort.10_pfaff_uncont fort.10_corr readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_agp_pfaff_uncont.dat
% cat corrsampling_agp_pfaff_uncont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468647E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.122986368E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.316227766E-07
Convert wavefunction from uncontracted JPf ansatz to contracted JPf ansatz:
cp fort.10_pfaff_uncont fort.10_in cp makefort10_pf_uncont.input makefort10_pf_cont.input
Edit makefort10_pf_cont.input
sed -i -e '/&symmetries/a nosym_contr=.true.' -e '/nshelljas/a ndet_hyb=2' makefort10_pf_cont.input
Convert wavefunction
makefort10.x < makefort10_pf_cont.input > out_make_pfaff_cont mv fort.10_new fort.10_out convertfort10-serial.x < convertfort10.input > out_conv_pfaff_cont cp fort.10_new fort.10_pfaff_cont
Clean wavefunction
cp fort.10_pfaff_cont fort.10 cleanfort10.x mv fort.10_clean fort.10 cp fort.10 fort.10_final cp fort.10 fort.10_corr
Run correlated sampling
cp fort.10_dft fort.10 readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_dft_pfaff_cont.dat
% cat corrsampling_dft_pfaff_cont.dat Number of bins 146 reference energy: E(fort.10) -0.148896680E+02 0.929468674E-02 reweighted energy: E(fort.10_corr) -0.148896681E+02 0.929468645E-02 reweighted difference: E(fort.10)-E(fort.10_corr) 0.115194734E-07 0.316227766E-07 Overlap square : (fort.10,fort.10_corr) 0.999999997E+00 0.505322962E-07
05-03 Li2 dimer and Li atom: VMC optimization (WF=JPf)
First, copy
fort.10andmakefort10.input.For Li2 dimer:
cp ../01convert_WF_JDFT_to_JAGP/fort.10_pfaff_cont_rot fort.10 cp fort.10 fort.10_sav cp ../01convert_WF_JDFT_to_JAGP/makefort10_pf_cont.input .
For Li atom:
cp ../01convert_WF_JDFT_to_JAGP/fort.10_pfaff_cont fort.10 cp fort.10 fort.10_sav cp ../01convert_WF_JDFT_to_JAGP/makefort10_pf_cont.input .
Run optimization
turborvb-serial.x < datasmin.input > out_min cleanfort10.x mv fort.10_clean fort.10
Check convergence and average
# check convergence readalles.x << ____EOS 1 1 0 1 0 1000 ____EOS # average readalles.x << ____EOS 1 81 1 0 0 1000 ____EOS
If you find many warnings in your output during optimization, do the following stabilization and continue optimization.
# Stabilization (especially needed for an open system, Li atom) cp fort.10 fort.10_in cp fort.10 fort.10_opt vi convertfort10mol.input # you should comment out epsdgm=-1d-14 -> !epsdgm=1d-14. # If this value is set negative, a generated molecular orbitals are random. # You can start from nmolmin=1, until nmolmin=N_el/ 2. # and check you do not loose much energy. convertfort10mol-serial.x < convertfort10mol.input > out_mol_stabilized cp fort.10_new fort.10
Run correlated sampling (aftre the above conversion)
cp fort.10_in fort.10_corr turborvb-serial.x < datasvmc.input > out_vmc readforward-serial.x < datasvmc.input > out_readforward mv corrsampling.dat corrsampling_pfaff_stabilized.dat cat corrsampling_pfaff_stabilized.dat
Convert back to a geminal (JAGP)
cp ../01convert_WF_JDFT_to_JAGP/makefort10_pf_cont.input . cp fort.10 fort.10_in set -- `grep -A1 "Parameters Jastrow two body" ./fort.10_in | tail -1` twobodyjas=$2 onebodyjas=$3 echo "twobodyjas=${twobodyjas} onebodyjas=${onebodyjas}" sed -i \ -e "/twobodypar/c\ twobodypar=${twobodyjas}" \ -e "/onebodypar/c\ onebodypar=${onebodyjas}" \ makefort10_pf_cont.input makefort10.x < makefort10_pf_cont.input > out_make_pfaff_cont mv fort.10_new fort.10_out convertfort10-serial.x < convertfort10.input > out_conv cp fort.10_new fort.10 cleanfort10.x mv fort.10_clean fort.10 cp fort.10_opt fort.10_new copyjas.x
Note
If you are using a more general spin-dependent Jastrow (i.e., twobody = -27), you should manually put Jastrow parameters after generating fort.10_out by makerfort10.x.
05-04 Li2 dimer and Li atom: VMC (WF=JPf)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
First, create a working directory, and copy the trial wavefunction:
cd ../03vmc/ cp ../02optimization/fort.10_new ./fort.10
Then, copy
datasmin.inputand rename it asave.in:cp ../02optimization/datasmin.input ave.in
You must also rewrite value of
ngenin ave.in asngen=1, andioptasiopt=1, by using an editor or typing the following commands:sed -i -e 's/ngen=.*/ngen=1/g' ave.in sed -i -e 's/iopt=.*/iopt=1/g' ave.in
Next, change
I/O flagin fort.10 to1which allow us to write the optimized variationa parameters:sed -i -e '/unconstrained/{n;s/0$/1/}' fort.10
Note
You may need to use GNU sed
gsedinstead ofsedon macOS.Run the dummy VMC by typing
turborvb-serial.x < ave.in > out_ave
Next step is to run VMC for calculating the total energy. Prepare
datasvmc.input, and run a VMC calculation by typing:turborvb-serial.x < datasvmc.input > out_vmc
After the VMC run finishes, check the total energy by running the script:
forcevmc.sh 10 2 1
A reblocked total energy and its variance is written in
pip0.d.
05-05 Li2 dimer and Li atom: LRDMC (WF=JPf)
Please refer to the Hydrogen tutorial for the details. Here, only needed commands are shown.
Prepare different working directories for each value of alat, copy fort.10 to each directory, and create an input file datasfn.input with the corresponding alat:
alat_0.8Z
alat_1.0Z
alat_1.2Z
alat_1.5Z
Then, run each LRDMC calculation after generating initial electron configurations at the VMC level.
cd ../04lrdmc/
cp ../03vmc/fort.10 .
lrdmc_root_dir=`pwd`
alat_list="0.8Z 1.0Z 1.2Z 1.5Z"
for alat in $alat_list
do
cd alat_${alat}
cp ${lrdmc_root_dir}/fort.10 ./fort.10
turborvb-serial.x < datasvmc.input > out_vmc;
turborvb-serial.x < datasfn.input > out_fn;
cd ${lrdmc_root_dir}
done
num=0
echo -n > ${lrdmc_root_dir}/evsa.gnu
for alat in $alat_list
do
cd alat_${alat}
num=`expr ${num} + 1`
echo "10 20 5 1" | readf.x
alat_d=`grep alat= datasfn.input | cut -f 2 -d '='`
echo -n "${alat_d} " >> ${lrdmc_root_dir}/evsa.gnu
tail -n 1 fort.20 | awk '{print $1, $2}' >> ${lrdmc_root_dir}/evsa.gnu
cd ${lrdmc_root_dir}
done
sed "1i 2 ${num} 4 1" evsa.gnu > evsa.in
funvsa.x < evsa.in > evsa.out
Finally we got:
Li2 dimer
HF = -14.8715 Ha
VMC(JDFT) = -14.9759(9) Ha
LRDMC(JDFT) = -14.989(1) Ha
VMC(JsAGPs) = -14.9821(4) Ha
LRDMC(JsAGPs) = -14.991(1) Ha
VMC(JAGPu) = -14.9819(5) Ha
LRDMC(JAGPu) = -14.990(1) Ha
VMC(JAGP) = -14.9831(4) Ha
LRDMC(JAGP) = -14.993(1) Ha
Exact = -14.9951 Ha
Li atom
HF = -7.4327 Ha
VMC(JDFT) = -7.4754(3) Ha
LRDMC(JDFT) = -7.4779(5) Ha
VMC(JsAGPs) = -7.477(2) Ha
LRDMC(JsAGPs) = -7.4783(4) Ha
VMC(JAGPu) = -7.477(3) Ha
LRDMC(JAGPu) = -7.4791(4) Ha
VMC(JAGP) = -7.4766(3) Ha
LRDMC(JAGP) = -7.4784(4) Ha
Exact = -7.4781 Ha
Binding energy
Ebond = 38.9 mHa (experiment)