# Tuning CTQMC Parameters for Nickel

# Issues with input and parameters

The aim of this tutorial is to give some indications on how to set the parameters of the DMFT calculation. We will give some examples of how the output should look like and, more importantly, how it shouldn’t.

The basic input files for the CTQMC solver are the *PARAMS* file, the *Trans.dat*, the *actqmc.cix* file, the *Delta.inp*, the *Eimp.inp* and finally the *status* files. Moreover **lmfdmft** requires the *indmfl* file. Descriptions of how to set the relevant parameters of the *PARAMS*, and the *indmfl* files are given below. You are not supposed to edit the other files, so we don’t give an explanation of them here.

Among the output files, the most important ones are *Sig.out* and *histogram.dat*. We will refer to these two to judge the quality of our calculation.

# Tuning input parameters for CTQMC+DMFT

Here we provide a standard *PARAMS* file for CT-QMC calculation for Nickel at =50 (eV) which should be run on at least 72 cores to produce high quality data over all energies.

```
OffDiagonal real
Sig Sig.out
Gf Gf.out
mode GH
svd_lmax 25
Delta Delta.inp
cix actqmc.cix
Nmax 3000 # Maximum perturbation order allowed
nom 200 # Number of Matsubara frequency points sampled
exe ctqmc # Name of the executable
tsample 50 # How often to record measurements
Ed [ -26.407760, -26.407519, -27.011761, -26.407397, -27.011682, -25.662041, -25.662114, -26.272431, -25.662012, -26.272482] # Ed: Eimp-Edc for PARAMS
M 80000000 # Total number of Monte Carlo steps per core
Ncout 2000000 # How often to print out info
PChangeOrder 0.9 # Ratio between trial steps: add-remove-a-kink / move-a-kink
mu 26.407760 # mu = -Eimp[0] for PARAMS
warmup 5000000 # Warmup number of QMC steps
sderiv 0.02 # Maximum derivative mismatch accepted for tail concatenation
aom 10 # Number of frequency points used to determin the value of sigma at nom
U 4.0
J 0.3
nf0 8.0
beta 50.0
```

The **mode GH** in the *PARAMS* file means that the Green’s functions (**G**) are sampled on Matsubara frequencies for first 200 frequency points as specified by the tage **nom 200** and the high energy tail is corrected by Hubbard solver (**H**). There are few other options available for the sampling modes like; **GM**, **SH**, **SM**. Here, **S** stands for self energy sampling and **M** stands for **analytic moment** corrections for the high energy tail.

Our experience says that GH is the most stable mode for a very good quality low energy and high energy data. However, these modes don’t solve the problem related to the noise at intermediate frequencies. To address that problem the expansion in SVD basis is introduced by Kristjan Haule. In his own language: **CTQMC was substantially reconfigured. The Green’s function or self-energy (or even the two particle vertex) can now be sampled both in the Matsubara frequency (old way) and in more efficient SVD basis ( These are SVD vectors of the analytic continuation kernel). This last method removes all noise at intermediate frequency.**

### A reasonably well converged dataset for Nickel for single-particle quantities:

Here you can find a very well converged dataset for Nickel for all the single-particle quantities from this link. The calculation converged after 14 iterations with a *PARAMS* file that look like above. To fasten the calculation what you can do is to add the **svd_lmax** switch only after 13 iterations, where you reliaze that all the quantities are converged, except for the fact that there is noise at intermediate frequencies. This calculation will run for 7-8 hours on 24 cores. You can, however, reduce the time by increasing the number cores and still acieve similar good statistics.

### Impose symmtery

However, if you want to impose symmetries within and sectors respectively, you want to modify the *indmfl.ni* file.

The *indmfl.ni* in that case looks like:

```
0.1 1.2 1 2 # hybridization Emin and Emax, measured from FS, renormalize for interstitials, projection type
1 0.0025 0.0025 600 -3.000000 1.000000 # matsubara, broadening-corr, broadening-noncorr, nomega, omega_min, omega_max (in eV)
1 # number of correlated atoms
1 1 0 # iatom, nL, locrot
2 2 1 # L, qsplit, cix
#================ # Siginds and crystal-field transformations for correlated orbitals ================
1 5 2 # Number of independent kcix blocks, max dimension, max num-independent-components
1 5 2 # cix-num, dimension, num-independent-components
#---------------- # Independent components are --------------
'x^2-y^2' 'z^2' 'xz' 'yz' 'xy'
#---------------- # Sigind follows --------------------------
1 0 0 0 0
0 1 0 0 0
0 0 2 0 0
0 0 0 1 0
0 0 0 0 2
#---------------- # Transformation matrix follows -----------
0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
0.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 -0.70710679 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000 0.00000000
-0.70710679 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70710679 0.00000000
```

### Switches for computing Two-particle Green’s functions

Till this point we have focussed only on the computation of single-particle quantities. However, if we want to compute spin and charge susceptibilities, the first thing we should have is the two-particle Green’s function. CT-QMC allows us to sample to two-particle quantities locally. To do so, we need to add two-switches to the *PARAMS* file.

```
SampleVertex 100
nOm 2
```

**SampleVertex 100** implies that a two-particle kink is added (removed) after every 100 kinks for the single-particle Green’s functions. In principle we are free to choose any integer for this particular switch, but 100 or more is recommended lest the computation time explodes. **nOm 2** stands for the number of bosonic frequencies on which we want to compute the two-particle Green’s functions. If someone is interested only in the lowest energy excitations for spin and charge, **nOm 1** would suffice. However, two-particle Green’s functions are explicit functions of two-fermionic frequencies and one bosonic frequency. If we don’t specify anything in the **PARAMS** file it is implied that the number of Fermionic frequencies **nomv** is 50.

### The *indmfl* file

*Note:* It is planned to change the syntax of this file. It will soon become obsolete!

This file is used by **lmfdmft** to define which atomic species and orbitals are mapped to the impurity. At the moment the file has a lot of redundant information, reminiscent of the formatting used in K. Haule’s DMFT package. Among all the information reported, you are interested only in a couple of variables at lines three and four.

```
1 # number of correlated atoms
1 1 0 # iatom, nL, locrot
```

The first variables defines how many atoms are correlated. *Note:* At the moment we are testing the code for values higher than 1. So far it hasn’t been possible to treat more than one atom.

The first variable of the second line (**iatom**=1 in the example) is the index of the correlated atom in the basis chosen. In the example, Ni has only one atom so it’s index is 1. This value has to be consistent with the *site.ext* file

If the site file is

```
ATOM=La POS= 0.5000000 0.5000000 -0.4807290
ATOM=La POS= -0.5000000 -0.5000000 0.4807290
ATOM=Cu POS= 0.0000000 0.0000000 0.0000000
ATOM=O POS= 0.0000000 0.5000000 0.0000000
ATOM=O POS= 0.0000000 0.0000000 0.6320273
ATOM=O POS= -0.5000000 0.0000000 0.0000000
ATOM=O POS= 0.0000000 0.0000000 -0.6320273
```

then **iatom=3**.

All the other values of the file are either ignored or should not be changed.

### The *PARAMS* file

This file is one of the input files read by the CTQMC sovler.

The variables contained in this file define the kind of calculation, allowing for a tuning of the Quantum Monte Carlo algorithm and details on how to treat the connection between the low-energy and the high-energy part of the self-energy. An example of the *PARAMS* file is reported in a dropdown box of the second tutorial.

**Basic parameters (U, J, nf0 and beta)**

Among the possible parameters are **U** and **J**, defining respectively the Hubbard in-site interaction and the Hund’s coupling constant in eV. **Warning:** The same value of **J** has to be passed to **atom_d.py** as well.

The variable **nf0** is the nominal occupancy of the correlated orbitals (e.g. **nf0 9** for , **nf0 8** for Ni).

Finally **beta** fixes the inverse temperature in eV.

**Warning:** Don’t forget to be consistent when you run **atom_d.py** (which wants **J** as an input) and **lmfdmft** (for which the flag **--ldadc=** has to be consistently defined as **U (nf0-0.5)-J(nf0-1)*0.5**).

**Setting the number of sampled frequencies (nom and nomD)**

The CTQMC solver gives a very accurate description of the self-energy in the low frequency range (for Matsubara’s frequencies close to 0), but it becomes too noisy at high frequencies.

Let be the total number of Matsubara’s frequencies. This number has been defined through **NOMEGA** in the *ctrl* file during the **lmfdmft** run. Only the first **nom** frequencies are actually sampled by the CTQMC solver, while the other points are obtained analytically through the approximated Hubbard 1 solver (high-frequency tail).

Excessively low values of **nom** will miss important features of the self-energy (e.g. a convex point in ), while too high values will give excessively noisy results. As a rule-of-thumb, a good initial guess is **nom** **beta**, but this value has usually to be adjusted during the DMFT loop. Some examples of how the looks like at different values of **nom** are given in the figures below.

**Setting the number of Monte Carlo steps (M and warmup)**

The higher is the number of Monte Carlo steps, the lower the noise in the QMC calculation. The parameter **M** defines the number of MC steps per core. Reliable calculations easily require at least 500 millions of steps in total. For instance, if you’re running on 10 cores, you can set **M 50000000**. You can judge the quality of your sampling by looking at the file *histogram.dat*. The closer it looks to a Gaussian distribution, the better is the sampling.

**Warning:** The variable **M** should be set keeping in mind that the higher it is, the longer the calculation. This is crucial when running on public clusters, where the elapsed time is computed per core. Too high values of **M** may consume your accounted hours very quickly! Moreover remember that you are supposed to broaden the output at each iteration, so you don’t actually need very clean *Sig.out*.

During the first **warmup** steps results are not accumulated, as it is normal on Monte Carlo procedures. This gives the ‘time’ to the algorithm to thermalise before the significative sampling. You can set **warmup**=**M**/1000.

**Setting the cutoff expansion order (Nmax)**

The variable **Nmax** defines the highest order accounted for in the hybrdization expansion. If you have chosen an excessively low values of **Nmax**, the *histogram.dat* file will be cut and will look weird, as shown below.

You should chose **Nmax** high enough for the Gaussian distribution of *histogram.dat* to be comfortably displayed. However note that the higher **Nmax** the longer the calculation, so chose values just above the higher Guassian tail.

**Warning:** the value of **beta** affects the number **Nmax**, so calculations on the same material at different temperatures will require different **Nmax**. At low **beta**, the Gaussian distribution is sharper and centered on lower order terms, as shown below. Therefore lower **beta** require lower **Nmax**.

**Connecting the tail (sderiv and aom)**

The connection between the QMC part and the Hubbard 1 part is done with a straight line starting at frequency number **nom**+1 and running until it intersect the Hubbard 1 self-energy. You can see it by comparing the *Sig.out* file with the *s_hb1.dat* (Hubbrad 1 only).

The two variables controlling the connection are **sderiv**, which is related to the angle of the straight line, and **aom** related to its starting point.

The straight line connecting the low- and the high-frequncy region can easily give rise to unwanted kinks. To broaden *Sig.out* is necessary to smooth these kinks out as well.

**Impurity levels (Ed and mu)**

The impurity levels reported at the fourth line of *Eimp.inp* enters in the *PARAMS* as the variable **Ed**. The variable **mu** is set as the first entry of the **Ed** with opposite sign.

**Warning:** These two variables change at every iteration so you have to constantly update their value throughout the DMFT loop.

### Phase transition boundaries

It may happen that, despite the high number of QMC steps, the *histogram.dat* file displays a double peak distribution simiar to the sum of two Gaussians. This is the case when the material is close to a phase transition.

In this case usually one or more channels of the self-energy are very noisy. One has to run for longer times, or use the status files to restart the calculation many times until only one peaks dominate and the histogram looks like a Gaussian. An example is given in the boxes below.

### Using status files

During the calculation, each core generates a *status* file. They contain some information about the sampling and should be used as restart files for other CTQMC calculations with similar parameters. They are read authomatically if they are in the folder where **ctqmc** is running.

They can be used basically in two ways.

- If you are performing iteration N, you should copy the
*status*files from iteration N-1 to speed up the convergence of the calculation. - If you realise that in one
**ctqmc**run, you haven’t achieved a good sampling (e.g.**M**too low, or close to phase transition), than you can run again the calculation without losing the effort already done. When restarted, the**ctqmc**solver will automatically read the*status*files and retrive their information.

**Warning:** Since there is one *status* file per core, you must pay attention to run on as many cores as *status* files you have. It should be safe (but not recommended) to run with a smaller number of cores, while running on more cores than *status* files gives wrong results.