• Lattice calculations
• Coordinate frame
• Converting lattice coordinate into array index
• Iterating over lattice points

Lattice calculations

W-SLDA Toolkit solves the problem on a Cartesian mesh grid with a lattice size N_x\times N_y\times N_z and lattice spacing D_x\times D_y\times D_z. The total extent of the simulation box is L_x\times L_y\times L_z=N_xD_x\times N_yD_y\times N_zD_z.
The lattice parameters must be provided at the compilation stage via predefines.h file:

/**
 * Define lattice size and lattice spacing
 * */
#define NX (8)
#define NY (10)
#define NZ (12)

#define DX (1.0)
#define DY (1.0)                                        
#define DZ (1.0)

In addition, we define parameters that you can use in all user-defined files:

// integration elements
#define DXYZ (DX*DY*DZ)
#define DXY (DX*DY)

// lattice points
#define NXYZ (NX*NY*NZ)
#define NXY (NX*NY)

// volume 
#define LX (DX*NX)
#define LY (DY*NY)
#define LZ (DZ*NZ)
#define LXYZ (LX*LY*LZ)
#define LXY (LX*LY)

Coordinate frame

Simulation box contains volume [0,DX*NX) x [0,DY*NY) x [0, DZ*NZ), where NX, NY, NZ and DX, DY, DZ are lattice dimensions and lattice spacings, respectively. The image below represents the simulation box for NX=NY=NZ=64 and DX=DY=DZ=1:

User-defined functions, like:

double v_ext(int ix, int iy, int iz, int it, int spin, double *params, size_t extra_data_size, void *extra_data)

take as input point coordinate (ix,iy,iz) expressed in lattice units. In order to convert it to Cartesian coordinates, use:

double x = DX*ix;
double y = DY*iy;
double z = DZ*iz;

Note: C language uses the first variable in the expression to set the precision of calculations, thus DX*ix will be executed using double precision (DX is double), while ix*DX will be executed using integer precision.
In many cases, it is convenient to set the domain origin to the center of the simulation box. This we do by shifting the frame by vector (NX/2, NY/2, NZ/2):

double x = DX*(ix-NX/2);
double y = DY*(iy-NY/2);
double z = DZ*(iz-NZ/2);

Converting lattice coordinate into array index

Variables (like density_a, delta, ...) are stored as one-dimensional arrays. Precisely, quantity for lattice coordinate (ix,iy,iz) is written in array element quantity[ixyz] where:

int ixyz = iz + NZ*iy + NZ*NY*ix;

Typically, access to quantities is provided via structures wslda_density and wslda_potential. These structures contain information about the coordinate frame that was used for storing the provided quantities:

typedef struct
{
    int nx;             
    int ny;             /// for 1d code it is set to 1
    int nz;             /// for 1d and 2d code it is set to 1
    int datadim;        /// dimensonality of data
    int blocklength;    /// number of elements in the array
...

It is recommended to use these variables, as they account for various dimensionality modes (3D, 2D, and 1D). For example, in order to extract the value of density in the center of the box, it is recommended to use:

// INPUT: wslda_density h_densities
int lNX=h_densities.nx, lNY=h_densities.ny, lNZ=h_densities.nz; // local sizes
int ix = lNX/2, iy = lNY/2, iz = lNZ/2; // (ix,iy,iz) points now to box center 
int ixyz = iz + lNZ*iy + lNZ*lNY*ix;
h_densities.density_a[ixyz]; // value of density in center of the box

Iterating over lattice points

The following scheme iterates over lattice points.

int lNX=h_densities.nx, lNY=h_densities.ny, lNZ=h_densities.nz; // local sizes
// or ...
// int lNX=NX, lNY=NY, lNZ=NZ;
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
    
    // ... compute something for coordinate (x,y,z) ....
    // rho_a[ixyz]=...;

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