Very slow for QMC sampling

[guuuj]

Hi developers,
I use cthyb impurity solver to do the DFT+DMFT calculations, and I found that for some systems, the calculation is very slow, even though I increase the temperature to beta=10.
My input is as follows:

[general]
seedname = "wannier90"
jobname = "dmft_dir"

beta = 10
n_iw = 401
n_tau = 4001
n_iter_dmft = 10

prec_mu = 0.002

g0_mix = 0.25

h_int_type = "kanamori"
h_int_basis = "qe"
U = [4.0, 4.0, 4.0]
J = [0.7, 0.7, 0.7]

dc_type = 0
dc_dmft = false

h5_save_freq = 5

calc_energies = true

mu_mix_const = 0.3
mu_mix_per_occupation_offset = 0.4
calc_mu_method = "brent"

[solver]
type = "cthyb"
length_cycle = 120
n_warmup_cycles = 4e+4
n_cycles_tot = 5e+6
imag_threshold = 1e-6

measure_density_matrix = true

n_l = 35
legendre_fit = true

[advanced]
dc_fixed_occ = [5, 5, 5]

I check the output, and I found that the acceptance of the QMC is very low. For example:

Warming up ...
11:46:05   0% ETA 00:29:21 cycle 2 of 40000
11:46:07   0% ETA 00:19:35 cycle 73 of 40000
11:46:10   0% ETA 00:18:26 cycle 169 of 40000
11:46:13   0% ETA 00:18:18 cycle 285 of 40000
11:46:17   1% ETA 00:18:43 cycle 418 of 40000
Insert error : recovering ... 
11:46:22   1% ETA 00:14:42 cycle 749 of 40000
11:46:28   2% ETA 00:13:35 cycle 1101 of 40000
11:46:36   3% ETA 00:13:01 cycle 1520 of 40000
11:46:46   4% ETA 00:13:07 cycle 1960 of 40000
Insert error : recovering ... 
11:46:58   6% ETA 00:13:16 cycle 2485 of 40000
11:47:13   7% ETA 00:13:34 cycle 3079 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
11:47:32  10% ETA 00:12:31 cycle 4150 of 40000
11:47:56  13% ETA 00:11:40 cycle 5462 of 40000
11:48:25  18% ETA 00:10:23 cycle 7353 of 40000
Insert error : recovering ... 
11:49:03  23% ETA 00:09:48 cycle 9268 of 40000
11:49:49  30% ETA 00:08:27 cycle 12241 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
11:50:47  38% ETA 00:07:25 cycle 15512 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
11:51:59  50% ETA 00:05:48 cycle 20180 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
11:53:30  65% ETA 00:03:54 cycle 26184 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
11:55:23  83% ETA 00:01:49 cycle 33461 of 40000
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
11:57:03 100% ETA 00:00:00 cycle 39999 of 40000
Insert error : recovering ... 



Accumulating ...
11:58:35   0% ETA 00:22:32 cycle 5 of 78125
Insert error : recovering ... 
11:58:37   0% ETA 00:23:31 cycle 117 of 78125
11:58:40   0% ETA 00:24:18 cycle 249 of 78125
11:58:43   0% ETA 00:21:57 cycle 462 of 78125
11:58:47   0% ETA 00:22:52 cycle 666 of 78125
11:58:52   1% ETA 00:23:31 cycle 918 of 78125
11:58:58   1% ETA 00:24:14 cycle 1216 of 78125
11:59:06   2% ETA 00:22:58 cycle 1705 of 78125
11:59:16   3% ETA 00:21:15 cycle 2401 of 78125
Insert error : recovering ... 
11:59:28   4% ETA 00:19:11 cycle 3412 of 78125
11:59:43   5% ETA 00:17:47 cycle 4662 of 78125
12:00:02   7% ETA 00:17:46 cycle 5870 of 78125
Insert error : recovering ... 
12:00:25   9% ETA 00:17:06 cycle 7586 of 78125
Insert error : recovering ... 
12:00:55  12% ETA 00:16:56 cycle 9456 of 78125
12:01:32  15% ETA 00:16:12 cycle 12024 of 78125
Insert error : recovering ... 
12:02:18  19% ETA 00:15:39 cycle 15000 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
12:03:16  24% ETA 00:14:45 cycle 18814 of 78125
12:04:28  30% ETA 00:13:38 cycle 23562 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
12:05:59  37% ETA 00:12:19 cycle 29287 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
12:07:52  46% ETA 00:10:44 cycle 36213 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
12:10:13  58% ETA 00:08:20 cycle 45505 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
12:13:09  73% ETA 00:05:14 cycle 57468 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
Insert error : recovering ... 
12:16:50  90% ETA 00:01:57 cycle 70556 of 78125
Insert error : recovering ... 
Insert error : recovering ... 
12:19:00 100% ETA 00:00:00 cycle 78124 of 78125
Insert error : recovering ... 


[Rank 0] Collect results: Waiting for all mpi-threads to finish accumulating...
[Rank 0] Timings for all measures:
Measure                                    | seconds   
Auto-correlation time                      | 0.168562  
Average order                              | 0.0373787 
Average sign                               | 0.0304049 
Density Matrix for local static observable | 38.2203   
G_tau measure                              | 4.45781   
Total measure time                         | 42.9144   
[Rank 0] Acceptance rate for all moves:
Move set Insert two operators: 0.00913899
  Move  Insert Delta_down_0: 0.00912923
  Move  Insert Delta_up_0: 0.00914875
Move set Remove two operators: 0.00913276
  Move  Remove Delta_down_0: 0.00911925
  Move  Remove Delta_up_0: 0.00914628
Move set Insert four operators: 0.000248019
  Move  Insert Delta_down_0_down_0: 0.000321701
  Move  Insert Delta_down_0_up_0: 0.000174047
  Move  Insert Delta_up_0_down_0: 0.000171026
  Move  Insert Delta_up_0_up_0: 0.000325317
Move set Remove four operators: 0.00025051
  Move  Remove Delta_down_0_down_0: 0.000327333
  Move  Remove Delta_down_0_up_0: 0.000175551
  Move  Remove Delta_up_0_down_0: 0.000173613
  Move  Remove Delta_up_0_up_0: 0.000325599
[Rank 0] Warmup lasted: 658.5 seconds [00:10:58]
[Rank 0] Simulation lasted: 1225.15 seconds [00:20:25]
[Rank 0] Number of measures: 78125
Total number of measures: 5000000
Average sign: 0.752246
Average order: 65.1337
Auto-correlation time: 249.875
Actual time for solver: 2049.41 s

I have set the rot_mat_type='hloc_diag' in Wannier90Converter, and I also test the parameter set_rot = 'hloc', but situation seems not change. I found that the calculation time is indeed system dependent, but the acceptance do not differ large for my another faster system.
So, my question is, is there any wrong in my settings? Is there any trick to accelerate the calculation? Or it is a normal case for difficult systems?

[alberto-carta]

Hi guuuj,

Could you provide more information on this system? How many orbitals is this?

Alberto

[the-hampel]

1. yes you should specify h_int_basis=‘wannier90' indeed. Note also that for this interaction U is not the intraorbital coulomb interaction term as in the Kanamori interaction, but rather the values of U=F0, J=(F2+F4)/14 .
2. actually n_cycles_tot controls directly the paramter n_cycles of cthyb and it is the number of measurements not sweeps. The name is a bit confusing. See here: https://triqs.github.io/cthyb/latest/guide/cthyb_convergence_tests.html for some extensive information of these parameters and also specifically https://triqs.github.io/cthyb/latest/guide/cthyb_convergence_tests.html#1.-length_cycle . So n_cycles_tot is the number of measurements, and length_cycle is the number of moves between two consecutive measurements.

[alberto-carta]

Hi guuj

is a 5 orbital impurity representing one of the TM atoms. And my system is CoO2, the Co4+ has 5 electrons. The calculations is a one-shot calculation. And my calculation is running on 64 cores. In fact, I learned a lot from your input, but I have some questions.

1. You mentioned that you have 4 TM atoms in your calculations, but I found that you specify 12 impurites in your input. I want to check that is it need to specify all the oxygen atoms in the input for a full d-p model? Or we can just focus on the TM atoms although we include oxygen orbitals in the wannier90_hr.dat?
2. You specify two keyword fit_min_n and fit_max_n. Are they important in the calculation? I don't know these keys before.
3. The most important thing is that you mention that we can restrict the local hilbert space by loc_n_min and loc_n_max, these are also two keywords that I don't know before. And it seems that it is important to accelerate the calculations. Are they also work during the warm up steps? It seems that in my case, the warm up is very hard. And how to determine the suitable values in my case? Does it restrict the nominal occupancy of the TM ion, or the actual occupancy?
4. For the hartree impurity, I notice that you set one_shot = true. In fact, I am not so understand this keyword. To compare the results with one-shot DMFT with cthyb solver, should we set one_shot = true?

I second @the-hampel that for a 5 orbital impurity for a pd model, U values are most often expressed in the Slater parametrization. For maximum efficiency I would suggest starting with h_int_type=density-density and move to full_slater later.
For best performance i would set rot_mat_type='hloc_diag in the wannier converter when constructing the input h5 file from w90, this will likely also make the sign problem better.

To answer your questions point by point

1. In the file that I sent there are two symmetry equivalent Cu atoms that are mapped to the same impurity and oxygen atoms that are treated explicitly at a simple hartree-fock level. I have a publication on this (https://www.nature.com/articles/s41524-025-01928-4), in general one would need to see case by case, but in many occasions a treatment of the ligand electrons is necessary to reproduce spectral and ground state features in pd models of materials.

2. I specify two keyword fit_min_n and fit_max_n. This is a relatively new feature that stabilizes the convergence of the calculation, it basically reduces the noise at high frequency of the green function. You can ignore it for now.

3. Yes they work also during warmup. I would look at the occupancy of your atom (say it has 6.5 electrons total in the shell), then i would set max excitations to 8 electrons (around two electrons more hop from the bath into the impurity) and the min to 4-5 electrons.
   Mind that the closer the number of electrons is to half filling, the more expensive the impurity problem becomes, since the number of available states at a given electron number grows as ( N_elec choose N_orbitals*2 )
   A higher or lower acceptance rate won't in general affect the runtime much, i think.

4. This is just not to spend too much time on converging the Hartree Fock impurity problem for the oxygens at the beginning of the calculation. You can safely ignore this flag for now in your calculations, but if at some point you start fiddling with the oxygens it might save you some computational time.

Anyway, in a pd model with 5 orbitals you can get away with only so much in terms of efficiency, this is a complicated problem.

Let me know if anything else is needed, for anything scientific you do not feel at ease sharing on Github feel free to contact me directly https://www.psi.ch/en/lms/people/alberto-carta

[guuuj]

Hi Alberto,
I really appreciate the offer to discuss scientific details directly. I’ll keep that in mind and reach out via email if I have specific questions that are better suited for private discussion. Here I have some further questions:
You said that the loc_n_min and loc_n_max can be determined by checking the occupancy of the atom. Is this occupancy the DFT one or DMFT one? At the first of dmft.out, there will print the onccupancy like Total density for imp 0 from DFT: , is this one? Or the one after some DMFT iterations summary of impurity observables: total occupany of impurity 0:. I found that in some cases, these values can differ by 1~2.
Best regards

[alberto-carta]

If you see that that is the case I would take into account both cases.

For instance if you Co goes from 4.5 electrons in DFT occupations (first total occupation in the observables) and at the end of the calculation this becomes 6 or 3 (last total occupation in the observables dat file), I would keep the loc_n_min to 2 and loc_n_max to 7. This is just an example

[guuuj]

I got it. Thank you for your quick reply. I will do some tests based on all your suggestions above.