• Imposing constraints on densities and potentials
  • Order of function calls
  • Example: imprinting vortex in the center of the system
• Constrained density functional theory calculations

Imposing constraints on densities and potentials

In each iteration, the user can impose by hand various constraints on densities or potentials. problem-definition.h file provides handlers for these operations.

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 arbitrarily modify 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
    }
}

Order of function calls

Functions are called as follows:

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 the 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 the 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
    }
}

Constrained density functional theory calculations

In the standard approach, we provide an external potential V_{ext}(\bm{r}), and after the minimization process, we obtain the density distribution. However, we can treat the external potential V_{ext}(\bm{r}) as a generalized Lagrange multiplier that controls the density. For example, we can ask: What is the shape of the external potential that provides the desired density distribution n_0(\bm{r})? Such a type of minimization is popular in nuclear physics and can also be implemented with API provided by W-SLDA Toolkit. In order to find such a potential, one can introduce a procedure similar to adjusting the chemical potential:

V_{ext}^{(it+1)}(\bm{r}) := V_{ext}^{(it)}(\bm{r}) + ab^{it}\frac{n^{(it)}(\bm{r}) - n_0(\bm{r})}{n_0(\bm{r}) + c}

where a,b,c are parameters that need to be adjusted to achieve convergence, it stands for the iteration number.

As an example, let us consider the quasi-1D case, where search for the external potential that produces density distribution in the form of two Gaussians:

n_0(\bm{r}) = A_1\exp\left(-\frac{(x-x_1)^2}{2\sigma_1^2}\right) + A_2\exp\left(-\frac{(x-x_2)^2}{2\sigma_2^2}\right)

In the figure below, the requested density is plotted by red, while with the blue color, we plot the density obtained after iterating the code with the prescription for the update of the external potential V_{ext}^{(it+1)}(\bm{r}):

The resulting external potential has a form as presented in the figure below: