Merge branch 'master' of https://github.com/computationalmodelling/fidimag

Author: marijanbeg

fidimag/atomistic/lib/util.c

@@ -8,10 +8,77 @@ inline double volume(double S[3], double Si[3], double Sj[3]) {
 	return tx + ty + tz;
 }
 
-// C = S_i \dot (S_{i+1} \times S_{j+1}) +  S_i \dot (S_{i-1} \times S_{j-1})
 double skyrmion_number(double *spin, double *charge,
                        int nx, int ny, int nz, int *ngbs) {
 
+    /* Calculation of the "Skyrmion number" Q for a two dimensional discrete
+     * spin lattice in the x-y plane (also known
+     * as finite spin chirality)
+     *
+     * The *spin array is the vector field for a two dimensional
+     * lattice with dimensions nx * ny
+     * (we can take a slice of a bulk from Python and pass it here, 
+     *  remember to do the ame for the neighbours matrix)
+     * The array follows the order:
+     *   [Sx0 Sy0 Sz0 Sx1 Sy1 Sz1 ... ]
+     *
+     * Charge is a scalar field array used to store the spin chirality /
+     * skyrmion number density (Q value per lattice site)
+     *
+     * *ngbs is the array with the neighbours information for every
+     * lattice site. The array is like:
+     *      [ 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, 1-x, 1+x, 1-y, ...  ]
+     *        i=0                           i=1                ...
+     *
+     * where  0-y  is the index of the neighbour of the 0th spin,
+     * in the -y direction, for example
+     *
+     *
+     * THEORY:
+     *
+     * Referring to the i-th site of a
+     * square lattice, we refer to the x direction for the i-th neighbours
+     * and to the y direction for the j-th neighbours, as here:
+     *
+     *                    X  (j + 1)
+     *                    |
+     *                (i) |
+     *     (i - 1) X------O------X   (i + 1)
+     *                    |
+     *                    |
+     *                    X   (j - 1)
+     *
+     * Then, we use the follwing expression:
+     *
+     * Q =  S_i \dot ( S_{i+1} \times S_{j+1} ) 
+     *      +  S_i \dot ( S_{i-1} \times S_{j-1} )
+     *
+     * This expression is based on the publication PRL 108, 017601 (2012)
+     * where Q is called "finite spin chirality". The idea comes from
+     * discrete chiral quantities in Hall effect studies. For example, at
+     * the end of page 3 in Rep. Prog. Phys. 78 (2015) 052502, it
+     * is argued:
+     *     scalar chirality (...) , which measures the volume enclosed 
+     *     by the three spins of the elementary triangle and, similarly to
+     *     (the vector chirlity) is sensitive to the sense of spin's 
+     *     rotation in the x–y plane
+     *
+     *  Hence we are taking the triangles formed by (i, i+1, j+1)
+     *  and (i, i-1, j-1) whose total area covers a unit cell,
+     *  and ommit the other two triangles (bottom right and top left)
+     *  When we sum this quantity across the whole lattice we cover every
+     *  triangle. The final quantity will be scaled by 8 PI
+     *  to match +-1 for a full skyrmion configuration
+     *
+     *  Recently, other ways to calculate a discrete skyrmion number have
+     *  been proposed: http://arxiv.org/pdf/1601.08212.pdf
+     *                 Phys. Rev. B 93, 024417 
+     *  
+     *  also based on using three spins using triangles. This could be
+     *  useful for applying to a hexagonal lattice in the future. 
+     *
+     */  
+
 	int i, j;
 	int index, id;
 
@@ -24,7 +91,9 @@ double skyrmion_number(double *spin, double *charge,
 	for (i = 0; i < nxy; i++) {
         index = 3 * i;
 
-        int idv = 6 * i;
+        /* The starting index of the nearest neighbours for the 
+         * i-th spin */
+        int id_nn = 6 * i;
 
         S[0] = spin[index];
         S[1] = spin[index + 1];
@@ -34,45 +103,52 @@ double skyrmion_number(double *spin, double *charge,
         S_j[0] = S_j[1] = S_j[2] = 0;
 
         // neighbour at -x
-        if (ngbs[idv] > 0) {
-            id = 3 * ngbs[idv];
+        // Remember that the index is -1 for sites without material
+        if (ngbs[id_nn] > 0) {
+            id = 3 * ngbs[id_nn];
             S_i[0] = spin[id];
             S_i[1] = spin[id + 1];
             S_i[2] = spin[id + 2];
         }
 
         // neighbour at -y
-        if (ngbs[idv + 2] > 0) {
-            id = 3 * ngbs[idv + 2];
+        if (ngbs[id_nn + 2] > 0) {
+            id = 3 * ngbs[id_nn + 2];
             S_j[0] = spin[id];
             S_j[1] = spin[id + 1];
             S_j[2] = spin[id + 2];
         }
 
+        // The  S_i \dot ( S_{i+1} \times S_{j+1} )
         charge[i] = volume(S, S_i, S_j);
 
         S_i[0] = S_i[1] = S_i[2] = 0;
         S_j[0] = S_j[1] = S_j[2] = 0;
 
         // neighbour at +x
-        if (ngbs[idv + 1] > 0) {
-            id = 3 * ngbs[idv + 1];
+        if (ngbs[id_nn + 1] > 0) {
+            id = 3 * ngbs[id_nn + 1];
             S_i[0] = spin[id];
             S_i[1] = spin[id + 1];
             S_i[2] = spin[id + 2];
         }
 
         // neighbour at +y
-        if (ngbs[idv + 3] > 0) {
-            id = 3 * ngbs[idv + 3];
+        if (ngbs[id_nn + 3] > 0) {
+            id = 3 * ngbs[id_nn + 3];
             S_j[0] = spin[id];
             S_j[1] = spin[id + 1];
             S_j[2] = spin[id + 2];
         }
 
+        //  The S_i \dot ( S_{i-1} \times S_{j-1} )
         charge[i]  += volume(S, S_i, S_j);
+
+        /* Scale the chirality quantity */
         charge[i] /= (8 * WIDE_PI);
 
+        /* We use the sum to output the total spin chirality
+         * or skyrmion number */
         sum += charge[i];
     }