In each iteration used can impose by hand various constraints on densities or potentials.

If you want to constrain density you need to provide a body of the function:

/**
 * THIS FUNCTION IS CALLED DURING THE SELF-CONSISTENT PROCESS.
 * Before each diagonalization process, user can modify arbitrarily densities
 * @param it iteration number
 * @param h_densities structure with densities, see (wiki) documentation for list of fields
 * @param params array of input parameters, before call of this routine the params array is processed by process_params() routine
 * @param extra_data_size size of extra_data in bytes, if extra_data size=0 the optional data is not uploaded
 * @param extra_data optional set of data uploaded by load_extra_data()
 * */
void modify_densities(int it, wslda_density h_densities, 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;
    
    // ITERATE OVER ALL POINTS
    ixyz=0;
    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
        
        // h_densities.rho_a[ixyz] stores value of spin-up particles densities for coordinate (x,y,z)
        // and similarly for other densities
        // ... below you can modify them at your wish ...
        
        
        ixyz++; // go to next point,  it should be last line of the triple loop
    }
}

If you want to constrain potentials you need to provide a body of the function:

/**
 * THIS FUNCTION IS CALLED DURING THE SELF-CONSISTENT PROCESS.
 * Before each diagonalization process, user can modify arbitrarily potentials
 * @param it iteration number
 * @param h_densities structure with densities, see (wiki) documentation for list of fields
 *                    NOTE: densities structure is processed by modify_densities(...) function before call ot this function.
 * @param h_potentials struture with potentials, see (wiki) documentation for list of fields
 * @param params array of input parameters, before call of this routine the params array is processed by process_params() routine
 * @param extra_data_size size of extra_data in bytes, if extra_data size=0 the optional data is not uploaded
 * @param extra_data optional set of data uploaded by load_extra_data()
 * */
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;
    
    // ITERATE OVER ALL POINTS
    ixyz=0;
    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
        
        // h_potentials.V_a[ixyz] stores value of spin-up particles mean-field potential for coordinate (x,y,z)
        // and similarly for other potentials
        // ... below you can modify them at your wish ...
        
        
        ixyz++; // go to next point,  it should be last line of the triple loop
    }
}

Ordering of function calls

Function are called as follow:

while( not (convergence criteria) )
{
    ...
    ... create a hamiltonian matrix (potentials as input) ...
    ... diagonalize the hamiltonian matrix ...
    ... compute densities ...
    modify_densities(it, densall, dc_params, extra_data_size, extra_data) ;
    ... compute potentials according to selected EDF ...
    modify_potentials(it, densall, potsall, dc_params, extra_data_size, extra_data) ;
    ....
}

Example: imprinting vortex in the center of the system

The standard method of imprinting quantum vortex is to impose that order parameter has the following structure:

\Delta(\vec{r})=|\Delta(\vec{r})|e^{i\phi}

where \phi is angle in cylindrical system, i.e \phi=\arctan(y/x). The procedure is implemented in function below:

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;
    
    // ITERATE OVER ALL POINTS
    ixyz=0;
    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
        
        // phase imprint, while absolute value of paring is adjusted automatically
        double d_abs = cabs(h_potentials.delta[ixyz]);  // take absolute value of paring field
        double phi = atan2(y,x);                      // compute by hand phase 
        h_potentials.delta[ixyz] = cexp(I*phi)*d_abs; // save new pattern of paring field

        ixyz++; // go to next point,  it should be last line of the triple loop
    }
}