First Tutorial: Learn how to use the CTQMC impurity solver

Details of the CTQMC solver

  
U  2.25             # Coulomb repulsion (F0)
mu 1.125            # Chemical potential
beta 100            # Inverse temperature
mode SM             # S stands for self-energy sampling, M stands for high frequency moment tail
M 5000000           # Number of Monte Carlo steps
svd_lmax 30         # number of SVD functions to project the solution  
nom 80              # number of Matsubara frequency points to sample
tsample 200         # how often to record the measurements
aom 1               # number of frequency points to determin high frequency tail
GlobalFlip 1000000  # how often to perform global flip
Delta Delta.inp     # Input bath function hybridization
cix one_band.imp    # Input file with atomic state

This is the file that enumerates all atomic/cluster states, in the absence of the fermionic bath, and also contains the matrix elements of the creation operators. Hence, one needs to perform an exact diagonalization of the atom to run the ctqmc. Note that this can be very useful, as the same code can run not just the Hubbard type models, but also Kondo type models, without any modification.

We first show the cix file for the single band model, which is already included in the iterate.py script:

# cluster_size, number of states, number of baths, maximum matrix size
1 4 2 1
# baths, dimension, symmetry, global flip
0       1 0 0
1       1 0 0
# cluster energies for unique baths, eps[k]
0 0
#   N   K   Sz size F^{+,dn}, F^{+,up}, Ea  S
1   0   0    0   1   2         3        0   0
2   1   0 -0.5   1   0         4        0   0.5
3   1   0  0.5   1   4         0        0   0.5
4   2   0    0   1   0         0        0   0
# matrix elements
1  2  1  1    1    # start-state,end-state, dim1, dim2, <2|F^{+,dn}|1>
1  3  1  1    1    # start-state,end-state, dim1, dim2, <3|F^{+,up}|1>
2  0  0  0
2  4  1  1   -1    # start-state,end-state, dim1, dim2, <4|F^{+,up}|2>
3  4  1  1    1
3  0  0  0
4  0  0  0
4  0  0  0
HB2                # Hubbard-I is used to determine high-frequency
# UCoulomb : (m1,s1) (m2,s2) (m3,s2) (m4,s1)  Uc[m1,m2,m3,m4]
0 0 0 0 0.0
# number of operators needed
0

Comment lines are ignored by the code, but the code expects some string or empty line. The second line contains the cluster size (1 for single site), the number of all atomic states (which is 4 for a single impurity, i.e., empty site, spin down, spin up, and doubly occupied site). This is followed by the number of bath blocks, i.e., the type of bands that hybridize with the impurity. Here we have spin-up and spin-down bath, hence number of baths is 2. Finally, the second line specifies the maximum size of the matrix, which defines the electron creation operator. In general, it is larger than one, when we have full Coulomb interaction.

Following the third line, we have several lines, which specify the type of baths we have. Since we specified that there are two baths, we need to have here two lines, one for spin-up and one for spin-down. The first number in the line is just increasing integer, the second specifies the dimension of the bath block. Each block is one dimensional here. Next we specify which blocks are equivalent due to symmetry. Here we have paramagnetic state, hence the two baths have the same index. Finally, we have a number, which specifies which baths should be flipped by the global flip. Two baths, which have the same index, will be exchanged.

After a comment line, we have single line for each atomic state (more precisely block of atomic states). Since we have four atomic states, there are four lines. First number stands for the index, second for site occupancy, third for momentum (which is irrelevant here), then it is Sz of the atomic state. This is followed by the dimension of the block. After that, we have so-called index table.

The idea of the cix file is that we block together all atomic eigenstates, which are coupled by creation operators (this is explained in PRB 75:155113 (2007) ). Once we have all blocks enumerates, we can have an index table, which says that creation operator for spin-down changes atomic block state 1 to state state 2, and creation operator for spin-up changes atomic state 1 to state 3. We also see that atomic state 2, under the action of creation operator for up-spin, results in doubly-occupied site (state 4), and that applying creation operator for down-spin is not allowed (Pauli principle). Hence, these two columns (in general there are as many columns as there are creation operators in the problem) specify the index table.

Finally, there are two more numbers sets in each line. The first stand for the atomic energy, and the second for the total spin quantum number. Here the energy excludes the trivial term due to Coulomb repulsion F0 and the chemical potential shift mu. Inside ctqmc code, each atomic state energy will be added a term U*N*(N-1)/2.-mu*N, where U and mu are specified in the params file, and N is specified in the second column.

Now that we specified the index table for creation operator, we need to specify the matrix elements of the creation operators. The following eight lines specify such matrix elements. In general there are number-of-atomic-state*number-of-baths lines.

For each atomic state, there are two lines, because we have two baths. The first is for spin-down and the second for spin-up creation operator. Here we need to follow the order chosen above in the index table. The first number in each line stands for the state on which creation operator acts, and the second is the resulting state number. These numbers should match the index table. The next two numbers specify the dimension of the starting and the resulting block (here everything is one dimensional). For resulting states, which are not allowed by the Pauli principle, we just write 0 for dimension. Finally, we give the matrix element value. In general we would have here a matrix of dim1*dim2, where dim1 is the dimension of the starting state, and dim2 the resulting state.

The next few lines specify the high frequency tails. In this version of the code, the key word HB1 or HB2 must follow the matrix elements. This forces ctqmc to either use Hubbard-1 (HB1) or computes moments analytically (HB2). But which one is chosen, can now be specified through mode=M or H in the PARAMS file, which overwrites this settings.

If we want to calculate moments analytically, i.e., mode=M, we need to specify the tensor of Coulomb repulsion. We need to specify here only the non-zero values of the Coulomb repulsion. Each line means one non-zero value for the four dimensional tensor of the Coulomb repulsion. However, we choose to exclude the trivial F0 term, which can be added later inside ctqmc code. Hence, for the single orbital model, there is no Coulomb U beyond F0, hence we do not need to specify anything. We choose to add a single line of zeros. In general, we would have four integer numbers for the index of for orbitals, and one floating point number for the value of the Coulomb repulsion.

Finally, we can sample any number of equal time operators during the run. The matrix elements of such operators would need to be specified here. We decided not to sample any other operator.

Since the single-orbital model is perhaps to trivial, we will explain also two other more complex examples. First, the two band model, and next a plaquette with four sites, and a single orbital per site. Let's start with the two band model, which can be downloaded here:

# Cix file for cluster DMFT with CTQMC
# cluster_size, number of states, number of baths
1 16 4 1
# baths, dimension, symmetry, global-flip
  0  1   0   0
  1  1   0   0
  2  1   1   1
  3  1   1   1
# cluster energies for unique baths, eps[k]
0 0 
#   N   K   Sz size
1   0   0    0  1     2   5   3   4      0   0   # 00
2   1   0  0.5  1     0   6   7   8      0   0   # u0
3   1   0  0.5  1     7  10   0   9      0   0   # 0u
4   1   0 -0.5  1     8  11   9   0      0   0   # 0d
5   1   0 -0.5  1     6   0  10  11      0   0   # d0
6   2   0    0  1     0   0  12  13      0   0   # 20
7   2   0    1  1     0  12   0  14      0   0   # uu
8   2   0    0  1     0  13  14   0      0   0   # ud
9   2   0    0  1    14  15   0   0      0   0   # 02
10  2   0    0  1    12   0   0  15      0   0   # du
11  2   0   -1  1    13   0  15   0      0   0   # dd
12  3   0  0.5  1     0   0   0  16      0   0   # 2u
13  3   0 -0.5  1     0   0  16   0      0   0   # 2d
14  3   0  0.5  1     0  16   0   0      0   0   # u2
15  3   0 -0.5  1    16   0   0   0      0   0   # d2
16  4   0    0  1     0   0   0   0      0   0   # 22
# matrix elements
1   2    1 1       1
1   5    1 1       1
1   3    1 1       1
1   4    1 1       1
2   0    0 0   
.....
.....
HB2
# UCoulomb : (m1,s1) (m2,s2) (m3,s2) (m4,s1)  Uc[m1,m2,m3,m4]
0  0  0  0    0
# number of operators needed
0

The four baths are all written in terms of a one dimensional blocks, and we assume that there is spin degeneracy, but the self-energy for the first orbital is different than for the second (symmetry index is different for the first two and the second two baths). Global flip exchanges the first two baths, and also the second two baths, since they have the same index.

Next we enumerate all states, and we give N, Sz, index table... for each state. After the comment, we also explain what each state stands for, for example up means the first orbital contains spin-up, and the second contains spin-down electron.

There are 16*4 lines of matrix elements, which we did not fully write out. One can examine them in the attached file.

Finally, let us give example of a plaquette (four-site) DMFT problem, with single orbital per-site. Here we need to perform the exact diagonalization with a chosen U=F0, because the cluster eigenvectors depend on the value of U, and the atomic state energies non-trivially depend on U. Hence, in the PARAMS file we will set U=0, and we will include U in the cix-file. We give here example for U=6t, which is close to the Mott transition.

  
# Cix file for cluster DMFT with CTQMC
# cluster_size, number of states, number of baths, maximum_matrix_size 
4 84 8 12
# baths, dimension, symmetry, global flip
0   1 0    0   # Sz=dn, K=(0,0)
1   1 0    0   # Sz=up, K=(0,0)
2   1 1    1   # Sz=dn, K=(pi,0)
3   1 1    1   # Sz=up, K=(pi,0)
4   1 1    2   # Sz=dn, K=(0,pi)
5   1 1    2   # Sz=up, K=(0,pi)
6   1 2    3   # Sz=dn, K=(pi,pi)
7   1 2    3   # Sz=up, K=(pi,pi)
# cluster energies for non-equivalent baths, eps[k]
-2 -0 2
 #   N   K   Sz size
 1   0   0    0   1      2   3   4   5   6   7   8   9                  0                  0 
 2   1   0 -0.5   1      0  10  11  12  14  15  17  18                 -2                0.5 
 3   1   0  0.5   1     10   0  12  13  15  16  18  19                 -2                0.5 
 4   1   1 -0.5   1     11  12   0  10  17  18  14  15                  0                0.5 
 5   1   1  0.5   1     12  13  10   0  18  19  15  16                  0                0.5 
 6   1   2 -0.5   1     14  15  17  18   0  10  11  12                  0                0.5 
 7   1   2  0.5   1     15  16  18  19  10   0  12  13                  0                0.5 
 8   1   3 -0.5   1     17  18  14  15  11  12   0  10                  2                0.5 
.......
# matrix elements
 1   2    1   1                     1 
 1   3    1   1                     1 
 1   4    1   1                     1 
 1   5    1   1                     1 
 1   6    1   1                     1 
 1   7    1   1                     1 
 1   8    1   1                     1 
 1   9    1   1                     1 
 2   0    0   0  
 2  10    1   4        0.952550727929                    0       0.188982236505      -0.238606003712 
 2  11    1   2                     1                    0 
......
HB1
# number of operators needed
0

Here we choose to work in the basis with momentum as a good quantum number, and we will choose baths to be enumerate by momentum and spin. We commented for each bath how was it chosen. Note that in the input/output file for ctqmc, we will find the same order of baths. But in Delta.inp we will specify only baths, which are different. Since we have only four different baths, we will give only for baths in Delta.inp file.

In the section on matrix elements, we have states of dimension higher than one. In this case, we need to specify dim1*dim2 numbers in a single row.

Finally, to perform the c-DMFT calculation for the Hubbard cluster, we give here an example of the python script iterate_cluster.py, which implements the DMFT self-consistency condition. To start, one also needs the cix-file hubbard_U_6_normal.cix. We also give the converged self-energy Sig.out, which should be the self-consistent solution at this U.

If one desires to calculate the cix file for other values of U, we want to give the tarbal with cix-file generator for the Hubbard model.

For good statistic, one should run the script on a parallel cluster with a hundred cores or so. One should specify in "mpi_prefix.dat" file, how should the parallel job be executed. For example "mpirun -n 100".

For completness we also mention more complicated example of d-wave superconducting state in plaquette DMFT.

# Cix file for cluster DMFT with CTQMC
# cluster_size, number of states, number of baths, maximum_matrix_size 
4 84 6 12
# baths, dimension, symmetry, global flip
  0   1       0                   0
  1   1     -0*                   1
  2   2       1    3    3  -1*    2
  3   2       1   -3   -3  -1*    3
  4   1       2                   4
  5   1     -2*                   5
# cluster energies for non-equivalent baths, eps[k]
-2 -0 2  0
 #   N   K   Sz size
 1   0   0    0   1      3   0   5   0   7   0   9   0                  0                  0 
 2   1   0 -0.5   1     10   1  12   0  15   0  18   0                 -2                0.5 
 3   1   0  0.5   1      0   0  13   0  16   0  19   0                 -2                0.5 
 4   1   1 -0.5   1     12   0  10   1  18   0  15   0                  0                0.5 
 5   1   1  0.5   1     13   0   0   0  19   0  16   0                  0                0.5 
 6   1   2 -0.5   1     15   0  18   0  10   1  12   0                  0                0.5 
 7   1   2  0.5   1     16   0  19   0   0   0  13   0                  0                0.5 
 8   1   3 -0.5   1     18   0  15   0  12   0  10   1                  2                0.5 

Useful output files

• nohup_imp.out.000
appears in parallel run (the same content in the serial run is printed to the screen), and can be monitored during the run to see the progress of the code.
• Sig.out
  The resulting impurity self-energy on Matsubara axis
• Gf.out
  The resulting impurity Green's function on Matsubara axis
• ctqmc.log
  Some log information on acceptance ratio and timings of various steps. It also contains
  • mom, which specifies the density matrix,
  • nf is the impurity occupancy,
  • chiS0 is the spin susceptibility at zero frequency,
  • chiD is the charge susceptibility at zero frequency,
  • <<Sz>> is the magnetization (in units of the spin, which needs to be multiplied by 2 to get units of bohr magneton).
  • Epot and Ekin stand for the potential and the kinetic energy. Notice that Epot contains not just the terms proportional to interaction U, but also the single particle terms due to crystal-field splittings. What we normally need, is the potential energy due to interaction only, which is named Tr(Sigma*G) or TrSigmaG and is mathematically equivalent to
    
  • Probability.dat
  The probability for all atomic states listed in cix file. The first column stands for the index of the block of atomic states (index of the superstate), the second column is the index within the block, the third column is the probability, which should be normalized to unity. Notice that with the probability of all atomic states, and information on atomic states, one can compute many physical properties. For example, the fluctuating moment <S^2> can be computed by taking S(S+1) from each atomic state, multiplied by its probability. Similarly, the potential energy can be obtained from atomic energies and their probabilities. The density matrix can be obtained by the density matrix of each atomic state, and their probabilities.
  • histogram.dat
  contains the probability for perturbation order k. Kinetic energy is <k> * T.
  • susc.dat
  Contains the spin and orbital susceptibility. But is printed only when "SampleSusc=1"
  • suscg.dat
  Contains the orbital susceptibility, but only when "SampleSusc=1" and the additional matrix elements are specified in cix file (see below).
  • Gcoeff.dat
  The coefficient in the SVD expansion of the Green's function in imaginary time. Relevant only when svd_lmax>0.
  • Gt.dat
  The imaginary time Green's function as expanded by the SVD functions. Relevant only when svd_lmax>0.
  • Gtau.dat
  Appears only when SampleGtau>0. The sparsely binned Green's function, which is of lower quality than Gt.dat or Gf.out, but contains raw data simply binned.
  • Gw.dat,Sw.dat
    The Fourier transform of Gt.dat and the self-energy (in mode=S in which self-energy is sampled directly). Relevant only when svd_lmax>0.

  The spin, the orbital and the charge susceptibility by CTQMC

  To calculate the z-component of spin-susceptibility
  
  or the charge susceptibility
  
  we use the matrix elements of Sz and N operator from cix file, and we evaluate the integral analytically for all times. To use this capability, we just need to add to PARAMS file the variable

  SampleSusc  1

  and the results appear in susc.dat. The first column is the spin susceptibility, and the second is the charge susceptibility.

  We can also control the number of sampled bosonic frequencies by nomb parameter.

  Notice that when susceptibility is sampled, we sample directly in Matsubara space and consequenctly switch svd_lmax=0.

  The orbital susceptibility is a bit more challenging to sample, as it requires to apply an operator, which is not necessary commuting with the local Hamiltonian. We implemented the formula, which samples any susceptibility in the subspace spanned by the partition function. Notice that if the operator requires jumping between blocks of the local Hamiltonian (such as for example x or y component of the spin-susceptibility), it is likely nonzero in the space outside the configuration space of the partition function, and then this implementation would give wrong answer. One would need to use the worm algorithm to sample a general susceptibility. If the susceptibility corresponding to an operator O is nonzero only in the configuration space of the partition function, we can compute

  
  where M^L contains all operators in the time evolution between t=0 and t=tau_i-1 and M^R contains operators between t=beta and t=tau_i. We store those time evolutions, and then we can quickly compute
  
  To use that feature, we need to specify the matrix elements of the O operator in the cix file. After the line

  # number of operators needed

  we normally specify a single number

  # number of operators needed
   1
  

  However, if we specify two numbers, like

  # number of operators needed
   1  1
  

  it means that we will first specify a local operator, for which we require an average value <O> and after that we will specify an operator O for which we want to compute susceptibility.

  The easiest way to use this feature is to run atom_d.py script, which performs exact diagonalization of the atom. We first prepare a file, which we call Trans.dat, that specifies the connection between the spherical harmonics and the orbitals used in ctqmc. Lets assume Trans.dat takes the form:

  5 5 #  size of Sigind and CF
  #---- Sigind follows
    1    0    0    0    0 
    0    2    0    0    0 
    0    0    3    0    0 
    0    0    0    4    0 
    0    0    0    0    5 
  #---- CF follows
    0.00000000   0.00000000    0.00000000   0.00000000    1.00000000   0.00000000    0.00000000   0.00000000    0.00000000   0.00000000 
    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.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.00000000  -0.70710679    0.00000000   0.00000000    0.00000000   0.00000000    0.00000000   0.00000000    0.00000000   0.70710679 
  

  In DFT+DMFT run, this information is normally stored in case.indmfl file, which you will learn in the following few tutorials. For now we just notice that Sigind is an index table for the form of the hybridization matrix, and CF is the unitary transformation from spherical to cubic harmonics. Each row stands for one cubic harmonics, and we use the order [z^2,x^2-y^2,xz,yz,xy] for the cubic harmonics, while the columns correspond to spherical harmonics with the order [-2,-1,0,1,2].

  Next we can execute

  atom_d.py J=0.8 l=2 OCA_G=False "CoulombF='Full'" "Eimp=[0.0,-0.1,-0.02,-0.02,-0.06]" "ORB=[0,0,1,-1,0,0,0,1,-1,0]"
  

  which means that we choose
  • Hund's coupling J=0.8eV (of the Slater form)
  • l=2 means d-orbitals (five orbital model)
  • OCA_G=False switches off the calculation of OCA input file
  • CoulombF='Full' specifies that we will use rotationally invariant Coulomb interaction
  • Eimp=[...] specifies the list of impurity levels
  • Finally, ORB specifies that this additional operator for orbital susceptibility should be printed in the cix file. The list contains 10 entries, as there are 10 orbitals+spins. Note that by convention we first specify 5 orbitals of spin-down, followed by the 5 orbitals of spin-up. The operator is constructed by
    
    For example, the specification ORB=[0,0,1,-1,0,0,0,1,-1,0] means
    
  Once we have such cix file with additional operator specified, we just run ctqmc with

  SampleSusc  1

  and the orbital susceptibility will appear in suscg.dat. The first column is the bosonic frequency, followed by the real and the imaginary part.

  Notes

  Finally, some further reading on details of the CTQMC implementation is available at these lecture notes.