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 a checkpoint generated by the code of a lower dimension
- Using a checkpoint generated on a lattice with a 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 already exists the user can read it from a file:

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

The 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 the values are then taken from the corresponding tags in the main part.

For spin-symmetric systems, the code typically converges without problems. For spin-imbalanced systems, the uniform solution is typically hard to obtain. To overcome this problem, we introduced a method that starts from the solution at a finite temperature init0Tstart and gradually decreases the temperature to the destination temperature init0Tstop with a 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 with default values, and the code indicates # WARNING: MAXITER REACHED!, indicating that self-consistency was not achieved. 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 also need to modify other init0 parameters.

NOTE: You can check the correctness of the settings of the uniform solver by using the 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 a checkpoint generated by the code of a lower dimension

The Toolkit allows initializing higher-dimensional calculations using checkpoints generated by lower-dimensional code; for example, 3D code can be initialized using 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 --- !!!
#        Dimensionality of the lattice has changed!
#        The code will change the dimensionality of the given checkpoint data to the new lattice.
# !!! --- ------- --- ------- --- ------- --- ------- --- ------- --- ------- --- !!!
# LOADING CHECKPOINT FILE `test/checkpoint.dat`
# CONVERTING CHECKPOINT FILE: 2D --> 3D

Using a checkpoint generated on a lattice with a different resolution.

If the provided checkpoint file is generated on a lattice with a 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 the value of potentials via the modify_potentials(...) routine. Initialization iteration has index it=-1. Below we present an 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 the 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 the 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 compromise the internal integrity of the data, which can affect the code's convergence properties. In such a case, it is suggested to disable mixing in the first iteration. When this option is applied, the data integrity will be retrieved after the first iteration.