Initialization of the st wslda solvers

The st-wslda codes support the following initialization modes:

inittype                0       #  0 - create uniform solution and start from it
                                # 10 - read uniform solution from file `inprefix`/uniform.solution
                                #  
                                #  5 - start from st-wslda checkpoint, inprefix points to folder with checkpoint binary file
                                #  
                                # -1 - custom initialization

Initializing solver with a uniform solution
Starting from checkpoint
- Using checkpoint generated by code of lower dimension
- Using checkpoint generated on a lattice with different resolution.
Custom initialization of the solver

Initializing solver with a uniform solution

The simplest way to start the self-consistent process is to initialize buffers with a solution of a uniform problem. To generate the uniform solution user should use:

inittype                0       #  0 - create uniform solution and start from it

Alternatively, if a uniform solution is already existing the user can read it from file:

inittype                10       #  10 - read uniform solution from file `inprefix`/uniform.solution 
inprefix              test       # uniform solution is stored in this location

Following flags in the input file are used to control the generation process of the uniform solution:

# ----------- IN CASE OF: inittype=0 ----------------
# init0Na                 28.0   # requested number of particles of type a, default=Na
# init0Nb                 28.0   # requested number of particles of type b, default=Nb
# init0muchange           0.1    # change rate of chemical potentials, default muchange
# init0Tstart             0.1    # start temperature, in units of eF, default value is equal to init0Tstop
# init0Tstop              0.01   # Temperature, in units of eF, default=temperature 
# init0DeltaT             0.01   # change of temperature in units of eF, default 0.01
# init0eps                1.0e-9 # epsilon for convergence, default 1.0e-9
# init0scmix              0.5    # mixing parameter in self-consistent process, default=linearmixing
# init0maxiter            1000   # maximum number of iterations, default=10000
init0debug              1      # debug level, default=0 (no debug info), 1 (basic debug info), 2 (detailed debug info)
init0save               1      # 0 - no saving of solution to file, 1 - save solution to the file

Note that, many of them are commented out, and then the values are taken from corresponding tags of the main code.

For spin symmetric systems to code typically converges without any problems. For spin-imbalanced systems, the uniform solution typically is hard to obtain. To overcome this problem, we introduced a method where we start from the solution at finite temperature init0Tstart and we gradually decrease the temperature to destination temperature init0Tstop with step init0DeltaT. Such a problematic case is shown below (for temperature 0.01):

# CREATING UNIFORM SOLUTION...
# DEBUG: n_a=0.017708, n_b=0.018750
# DEBUG: eF_a=0.516086, eF_b=0.536132, eF_avg=0.526157
# DEBUG: N_a=17.000000, N_b=18.000000
# DEBUG: p=-0.028571
# DEBUG: alph_a=1.089111, alph_b=1.098021, alph_plus=1.093566
# DEBUG: dalphm_dna=4.393928, dalphm_dnb=-4.149821
# DEBUG: dalphp_dna=0.856253, dalphp_dnb=-0.808683
# DEBUG: dtildeC_dna=-0.323965, dtildeC_dnb=-0.274273
# DEBUG: dD_dna=-0.289825, dD_dnb=-0.265175
# DEBUG: D=-0.006063, tC=-0.032638
# DEBUG: tau_a=0.010967, tau_b=0.012063
# DEBUG: delta=0.258043, nu=0.007702
# TEMPCONV: T=0.010000, iter=10000, delta/eF_a=0.509012, mu_a/eF_a=-0.061340, delta/eF_b=0.489980, mu_b/eF_b=0.911447
# WARNING: MAXITER REACHED!
# TEMPCONV: T=0.010000, energy_kin=1.923270, energy_pot=-0.526501, energy_pair=-0.967234, energy_tot=0.429535
# ENERGY CUT-OFF: EC1=5.150865
# ENERGY CUT-OFF: EC2=-4.682460
# ENERGY CUT-OFF: EC=4.909640
# UNIFORM SOLUTION: delta/eF_a=  0.5090, mu_a/eF_a= -0.0613, delta/eF_b=  0.4900, mu_b/eF_b=  0.9114, ec=  4.9096
# UNIFORM SOLUTION: energy_kin=  1.923270037681, energy_pot= -0.526501268447, energy_pair= -0.967234017120, energy_tot=  0.429534752115
# UNIFORM SOLUTION: nwf=970

It was obtained for default values, and the code indicates # WARNING: MAXITER REACHED! which means that the self-consistency was not obtained. However, when setting:

init0Tstart             0.1    # start temperature, in units of eF,

we obtain:

# CREATING UNIFORM SOLUTION...
# DEBUG: n_a=0.017708, n_b=0.018750
# DEBUG: eF_a=0.516086, eF_b=0.536132, eF_avg=0.526157
# DEBUG: N_a=17.000000, N_b=18.000000
# DEBUG: p=-0.028571
# DEBUG: alph_a=1.089111, alph_b=1.098021, alph_plus=1.093566
# DEBUG: dalphm_dna=4.393928, dalphm_dnb=-4.149821
# DEBUG: dalphp_dna=0.856253, dalphp_dnb=-0.808683
# DEBUG: dtildeC_dna=-0.323965, dtildeC_dnb=-0.274273
# DEBUG: dD_dna=-0.289825, dD_dnb=-0.265175
# DEBUG: D=-0.006063, tC=-0.032638
# DEBUG: tau_a=0.010967, tau_b=0.012063
# DEBUG: delta=0.258043, nu=0.007702
# TEMPCONV: T=0.100000, iter=2911, delta/eF_a=0.503474, mu_a/eF_a=0.195268, delta/eF_b=0.484650, mu_b/eF_b=0.668692
# TEMPCONV: T=0.100000, energy_kin=1.915994, energy_pot=-0.526501, energy_pair=-0.948608, energy_tot=0.440885
# TEMPCONV: T=0.090000, iter=1873, delta/eF_a=0.503667, mu_a/eF_a=0.163986, delta/eF_b=0.484835, mu_b/eF_b=0.698980
# TEMPCONV: T=0.090000, energy_kin=1.916098, energy_pot=-0.526501, energy_pair=-0.949330, energy_tot=0.440267
# TEMPCONV: T=0.080000, iter=1482, delta/eF_a=0.503698, mu_a/eF_a=0.132936, delta/eF_b=0.484865, mu_b/eF_b=0.729129
# TEMPCONV: T=0.080000, energy_kin=1.915824, energy_pot=-0.526501, energy_pair=-0.949435, energy_tot=0.439888
# TEMPCONV: T=0.070000, iter=1130, delta/eF_a=0.503659, mu_a/eF_a=0.102453, delta/eF_b=0.484828, mu_b/eF_b=0.758759
# TEMPCONV: T=0.070000, energy_kin=1.915372, energy_pot=-0.526501, energy_pair=-0.949276, energy_tot=0.439595
# TEMPCONV: T=0.060000, iter=822, delta/eF_a=0.503586, mu_a/eF_a=0.072729, delta/eF_b=0.484758, mu_b/eF_b=0.787641
# TEMPCONV: T=0.060000, energy_kin=1.914814, energy_pot=-0.526501, energy_pair=-0.948989, energy_tot=0.439323
# TEMPCONV: T=0.050000, iter=561, delta/eF_a=0.503481, mu_a/eF_a=0.043927, delta/eF_b=0.484656, mu_b/eF_b=0.815579
# TEMPCONV: T=0.050000, energy_kin=1.914134, energy_pot=-0.526501, energy_pair=-0.948580, energy_tot=0.439052
# TEMPCONV: T=0.040000, iter=351, delta/eF_a=0.503324, mu_a/eF_a=0.016246, delta/eF_b=0.484505, mu_b/eF_b=0.842323
# TEMPCONV: T=0.040000, energy_kin=1.913266, energy_pot=-0.526501, energy_pair=-0.947983, energy_tot=0.438782
# TEMPCONV: T=0.030000, iter=223, delta/eF_a=0.503090, mu_a/eF_a=-0.010016, delta/eF_b=0.484280, mu_b/eF_b=0.867514
# TEMPCONV: T=0.030000, energy_kin=1.912119, energy_pot=-0.526501, energy_pair=-0.947102, energy_tot=0.438515
# TEMPCONV: T=0.020000, iter=221, delta/eF_a=0.502800, mu_a/eF_a=-0.034271, delta/eF_b=0.484001, mu_b/eF_b=0.890670
# TEMPCONV: T=0.020000, energy_kin=1.910759, energy_pot=-0.526501, energy_pair=-0.946017, energy_tot=0.438242
# TEMPCONV: T=0.010000, iter=112, delta/eF_a=0.502792, mu_a/eF_a=-0.054949, delta/eF_b=0.483994, mu_b/eF_b=0.911258
# TEMPCONV: T=0.010000, energy_kin=1.910379, energy_pot=-0.526501, energy_pair=-0.945957, energy_tot=0.437921
# ENERGY CUT-OFF: EC1=5.146473
# ENERGY CUT-OFF: EC2=-4.682285
# ENERGY CUT-OFF: EC=4.907524
# UNIFORM SOLUTION: delta/eF_a=  0.5028, mu_a/eF_a= -0.0549, delta/eF_b=  0.4840, mu_b/eF_b=  0.9113, ec=  4.9075
# UNIFORM SOLUTION: energy_kin=  1.910378825876, energy_pot= -0.526501268447, energy_pair= -0.945956681969, energy_tot=  0.437920875460
# UNIFORM SOLUTION: nwf=970

Depending on the settings you will need to also modify other init0 parameters.

NOTE: You can check the correctness of settings of the uniform solver by using solve-uniform tool. For more info see here.

Starting from checkpoint

The code can be initialized by an already existing solution generated by st-wslda. It is the most commonly used option. To activate this option you need to set:

inittype                5       #  5 - start from st-wslda checkpoint, inprefix points to folder with checkpoint binary file
inprefix             test

Using checkpoint generated by code of lower dimension

The Toolkit allows for initializing higher dimensional calculations via checkpoint generated by lower-dimensional code, for example, 3D code can be initialized by checkpoints generated by 2D or 1D code.
In the case of resizing the checkpoint file, a message of this type will be provided on stdout:

# CODE: ST-WSLDA-3D
# VERSION: 2021.02.22
# LATTICE: 8 x 10 x 12
...
# INSPECTING CHECKPOINT FILE `test/checkpoint.dat`
# CHECKPOINT FOR 2D LATTICE: [NX,NY,NZ]=[8,10,12], [DX,DY,DZ]=[1.000,1.000,1.000], [LX,LY,LZ]=[8.000,10.000,12.000]
# !!! --- WARNING --- WARNING --- WARNING --- WARNING --- WARNING --- WARNING --- !!!
#        Dimensonality of the lattice has changed!
#        The code will change dimensionality of given checkpoint data to the new lattice.
# !!! --- ------- --- ------- --- ------- --- ------- --- ------- --- ------- --- !!!
# LOADING CHECKPOINT FILE `test/checkpoint.dat`
# CONVERTING CHECKPOINT FILE: 2D --> 3D

Using checkpoint generated on a lattice with different resolution.

If provided checkpoint file is generated on a lattice with different resolution, than specified in predefines.h, then automatic interpolation will be applied.

Custom initialization of the solver

The static codes can start the self-consistent process from arbitrary initial conditions. To activate this option user must specify in the input file:

inittype                -1       # -1 - custom initialization

Using this functionality it is mandatory to initialize value of potentials via modify_potentials(...) routine. Initializtion iteration has index it=-1. Below we present example of such initialization.

void modify_potentials(int it, wslda_density h_densities, wslda_potential h_potentials, double *params, size_t extra_data_size, void *extra_data)
{    
    // DETERMINE LOCAL SIZES OF ARRAYS (CODE DIMENSIONALITY DEPENDENT)
    int lNX=h_densities.nx, lNY=h_densities.ny, lNZ=h_densities.nz; // local sizes
    int ix, iy, iz, ixyz;
    
    if(it==-1 && wsldapid==0) wprintf("SETTING MY VALUES FOR STARTING POINT\n");
    
    // ITERATE OVER ALL POINTS
    ixyz=0;
    if(it==-1) for(ix=0; ix<lNX; ix++) for(iy=0; iy<lNY; iy++) for(iz=0; iz<lNZ; iz++)
    {
        double x = DX*(ix-lNX/2);
        double y = DY*(iy-lNY/2); // for 1d code y will be always 0
        double z = DZ*(iz-lNZ/2); // for 1d and 2d codes z will be always 0
        
        // my custom initialization of the solver 
        
        h_densities.rho_a[ixyz]=h_densities.rho_b[ixyz]=0.018; // guess for the density, here I set constant value 
        
        // all potentials must be filled, since they are used to construct BdG matrix
        h_potentials.delta[ixyz]=0.1 + I*0.0; // guess for delta
        h_potentials.alpha_a[ixyz]=h_potentials.alpha_b[ixyz]=1.0; // bare mass
        h_potentials.V_a[ixyz]=h_potentials.V_b[ixyz]=0.0; // no mean-field
        h_potentials.A_a_x[ixyz]=h_potentials.A_b_x[ixyz]=0.0; // no current potentials
        h_potentials.A_a_y[ixyz]=h_potentials.A_b_y[ixyz]=0.0; // no current potentials
        h_potentials.A_a_z[ixyz]=h_potentials.A_b_z[ixyz]=0.0; // no current potentials
        
        ixyz++; // go to next point,  it should be last line of the triple loop
    }
    
    // my guess for chemical potentials 
    if(it==-1)
    {
        h_potentials.mu[SPINA]=0.5;
        h_potentials.mu[SPINB]=0.5;
    }
}

Note, that when you initialize the solver manually you may break the internal integrity of data, and it can impact the convergence properties of the code. In such case, it is suggested to disable mixing in the first iteration. When applied this option, the integrity of data will be retrieved after the first iteration.