Convergence control

Convergence criteria

There are two parameters controlling the computation accuracy of the solution (listed in input file):

energyconveps           1.0e-6   # energy convergence epsilon
                                 # energy changes between two subsequent iterations 
                                 # expressed in units, E_ffg needs to be smaller than npartconveps
                                 # for converged solution
                                 # default=1.0e-6
npartconveps            1.0e-6   # number of particles convergence epsilon
                                 # relative difference between obtained particle number
                                 # and total particle number must be smaller than npartconveps
                                 # for converged solution
                                 # default=1.0e-6
# npartconveps_a          1.0e-6   # or you can control npartconveps for each component separately using tags with `a` ...
# npartconveps_b          1.0e-6   # ... and `b` suffixes.

The solution is regarded as converged when all contributions to the energy satisfy energyconvep for each entry (E_{\textrm{ffg}}=\frac{3}{5}N\varepsilon_F):

\dfrac{|E_{\textrm{entry}}^{(i)}-E_{\textrm{entry}}^{(i-1)}|}{E_{\textrm{ffg}}}<\varepsilon_{\textrm{energy}}

and particle number satisfies criteria npartconveps:

\dfrac{|N_{\sigma}^{(i)}-N_{\sigma}^{(input)}|}{N_{\uparrow}^{(input)}+N_{\downarrow}^{(input)}}<\varepsilon_{\textrm{npart}}^{(\sigma)}

Info about the convergence status is printed to standard output. Example of a report for converged solution for input parameters:

Na                      17       # requested particle number, spin-a
Nb                      18       # requested particle number, spin-b
energyconveps           1.0e-6   # energy convergence epsilon
npartconveps            1.0e-6   # number of particles convergence epsilon

may look like:

# CONVERGENCE REPORT PARTICLE NUMBER: it=58
    SPINA: NEW=              17 OLD=       17.000003 DIFF=  -3.1644056e-06 CONVSTATUS=  PASS NPARTCONV=   2.1975596e-09
    SPINB: NEW=       18.000001 OLD=       17.999995 DIFF=   5.4727832e-06 CONVSTATUS=  PASS NPARTCONV=    2.406396e-08
# CONVERGENCE REPORT ENERGY: it=58
    E_kin: NEW=        1.879979 OLD=       1.8799781 DIFF=   8.9986865e-07 CONVSTATUS=  PASS
    E_pot: NEW=     -0.51360103 OLD=     -0.51360098 DIFF=  -4.7134213e-08 CONVSTATUS=  PASS
   E_pair: NEW=     -0.94237189 OLD=     -0.94237177 DIFF=  -1.2408755e-07 CONVSTATUS=  PASS
   E_curr: NEW=   3.5892537e-14 OLD=   9.6882531e-14 DIFF=  -6.0989994e-14 CONVSTATUS=  PASS
 E_potext: NEW=               0 OLD=               0 DIFF=               0 CONVSTATUS=  PASS
E_pairext: NEW=               0 OLD=               0 DIFF=               0 CONVSTATUS=  PASS
 E_velext: NEW=               0 OLD=               0 DIFF=               0 CONVSTATUS=  PASS
  ------------------------------------------------------------------------------------------
    E_tot: NEW=      0.42400609 OLD=      0.42400536 DIFF=   7.2864682e-07 CONVSTATUS=  PASS

Number of iterations

The algorithm iterates until one of the following criteria is met:

Convergence is achieved (see the section above).
Number of iterations exceeded the maximum allowed number of iterations specified in the input file

maxiters                200     # maximum number of iterations, default=10000

Imposing constraints to speed up convergence

Providing additional information about the expected solution's properties (based on physical arguments) can significantly speed up convergence.

spinsymmetry            1       # set to 1 if you want to impose Na=Nb, default: 0
nocurrents              1       # set to 1 if you want to impose ja=jb=0 
                                # (solution has no currents), default: 0

Mixing types

The self-consistent algorithm tracks the behavior of solutions as a function of iteration number and uses this information to generate a new prediction for the solution. In the samples scheme (called linear mixing), this process is realized via the formula:

for(ixyz = 0; ixyz < SOLDIM; ixyz++) 
  h_solution[ixyz] = md.linearmixing * h_solution[ixyz] + (1.0-md.linearmixing) * h_solution_old[ixyz];

Here, there is a subtle issue: what should be treated as the solution vector: densities of potentials? In principle, this should not matter, but in practice, it turns out that convergent properties depend on this choice. For this reason W-SLDA toolkit provided a tag in the input file that can be used to selectthe meaning of the solution vector:

mixingtype              p   # 'd' - mix densities, 'p' - mix potentials (default)

Our tests show that the algorithm typically converges faster when potentials are used in the mixing procedure. However, we found (rare) cases where mixing densities improves efficiency.

Disabling mixing in the first iteration

In some cases, it is required to disable mixing in the first iteration. Typically, it is required when:

The solver has been initialized manually, and thus the integrity of data is not maintained.
automatic interpolations have been applied.

In such cases, not all buffers related to the previous iteration are filled correctly. In the case of the custom initialization, most likely they are set to zero, while in the case of the automatic interpolations, some quantities like \nu or \tau are interpolated incorrectly (strong cut-off dependence is not taken into account). To prevent incomplete or incorrect data from propagating into the self-consistent process in the first iteration, we can disable code fragments that depend on historical iterations. To do this, use the following flag:

nomixstart              1       # do not do mixing for starting iteration