Constraining-densities-and-potentials.md

+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:
+```c
+/**
+ * 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:
+```c
+/**
+ * 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:
+```c
+while( not (convergence criteria) )
+{
+    ... create hamiltoin matrix (potententials as input) ...
+
+    ... diagonalize hamitonian ...
+    
+    ... 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:  
+```math
+\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:
+
+
+```c
+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
+    }
+}
+```
\ No newline at end of file