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 iteration 
                                 # 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 status of convergence is printed on 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 iteration exceeded maximal allowed number of iteration specified in input file

maxiters                200     # maximum number of iterations, default=10000

Imposing constraints to speed up convergence

Providing additional information about properties of the expected solution (based on physical arguments) may speed-up converging significantly.

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 if the iteration number and utilizes this information to create 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 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 converges properties depend on this choice. For this reason W-SLDA toolkit provided a tag in the input file that can be used to select meaning of the solution vector:

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

Our tests show that typically the algorithm converges faster when potentials are used for the mixing procedure. However, we found (rare) cases where mixing of densities provides better 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 integral 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 do not propagate the incomplete/incorrect data into the self-consistent process in the first iteration we can disable fragments of the code that depend on the historical iterations. To do this use the following flag:

nomixstart              1       # do not do mixing for starting iteration