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 extend 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 mast be provided at 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. 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 coordinate use:

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

Note: C language uses first variable in expression to set 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 origin of simulation of domain in center of the simulations 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 provided there 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 quantity array = nx*ny*nz
...

It is recommended to use these variables, as they take into account various dimensionality modes (3D, 2D and 1D). For example, in order to extract value of density in 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 is used to iterate 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
}