• 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 the body of the function:

/**
 * THIS FUNCTION IS CALLED DURING THE SELF-CONSISTENT PROCESS.
 * Before each diagonalization process, the user can modify densities arbitrarily.
 * @param it iteration number
 * @param h_densities structure with densities, see (wiki) documentation for the list of fields
 * @param params array of input parameters, before the 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 the next point, it should be the 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, the user can modify potentials arbitrarily.
 * @param it iteration number
 * @param h_densities structure with densities, see (wiki) documentation for the list of fields
 *                    NOTE: densities structure is processed by modify_densities(...) function before call of this function.
 * @param h_potentials structure with potentials, see (wiki) documentation for list of fields
 * @param params array of input parameters, before the 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 the next point, it should be the 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) ;
    ....
}

For the complete workflow, see here.

Example: imprinting vortex in the center of the system

The standard method of imprinting a 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 the 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 the API provided by the 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 we search for the external potential that produces a 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 in red, while the blue color we to 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: Inspect the attached files to learn about the protocol's example implementation. The most important functions used in the implementation are

• get_extra_data_size(...): prepares buffers for the external potential and the desired shape of the density.
• load_extra_data(...): initializes the guess for the external potential and shape of the density.
• modify_potentials(...): implements the external potential update protocol.
• v_ext(...): returns the value of the external potential in the given iteration.

Files:

• rhoconstr_predefines.h
• rhoconstr_problem-definition.h
• input.txt
• rhoconstr.stdout