Documented and added methods for the calculation of discrete skyrmion number or topological charges. The LLG code was modified to check which meshes are appropriate for every method

Author: davidcorteso

fidimag/atomistic/lib/util.c

@@ -159,6 +159,11 @@ double skyrmion_number(double *spin, double *charge,
 }
 
 double dot(double *a, double *b){
+    /* Dot product for vectors of 3 components, given by arrays a and b.  If
+     * pointers are passed, it uses the first three components starting from
+     * that place in memory, i.e. if we pass &a[2], where a = {0, 1, 2, 3, 4},
+     * the dot product uses {2, 3, 4} as a vector (same for b)
+     */
     double dp = 0;
     int i;
 
@@ -167,6 +172,35 @@ double dot(double *a, double *b){
 }
 
 double compute_BergLuscher_angle(double *s1, double *s2, double *s3){
+
+    /* Compute the spherical angle given by the neighbouring spin 3-vectors s1,
+     * s2 and s3 (these are arrays that can be also passed as pointers; see the
+     * -dot- function). The spherical angle Omega, defined by the
+     *  vectors in a unit sphere, is computed using the imaginary exponential
+     *  defined by Berg and Luscher [Nucl Phys B 190, 412 (1981)]:
+
+     *   exp (i sigma Omega) = 1 + s1 * s2 + s2 * s3 + s3 * s1 + i s1 * (s2 X s3)
+     *                         ------------------------------------------------
+     *                            2 (1 + s1 * s2) (1 + s2 * s3) (1 + s3 * s1)
+
+     * The denominator is a normalisation factor which we call rho, and i
+     * stands for an imaginary number in the numerator. The factor sigma is an
+     * orientation given by sign(s1 * s2 X s3), however, we do not use it since
+     * we take counter clock wise directions for the (s1, s2, s3) triangle of
+     * spins, as pointed out by Yin et al. [PRB 93, 174403 (2016)].
+     *
+     * Therefore we use a complex logarithm:
+
+     *      clog( r * exp(i theta) ) = r + i theta
+
+     * to calculate the angle Omega, since the clog is well defined in the
+     * [-PI, PI] range, giving the correct sign for the topological number (we
+     * could also use the arcsin when decomposing the exp).
+     *
+     * Notice we normalise the angle by a 4 PI factor
+     *
+     */
+
     double rho;
     double complex exp;
     double crossp[3];
@@ -180,8 +214,8 @@ double compute_BergLuscher_angle(double *s1, double *s2, double *s3){
                  * (1 + dot(&s3[0], &s1[0]))
                );
 
-    exp = (1 + dot(&s1[0], &s2[0]) 
-             + dot(&s2[0], &s3[0]) 
+    exp = (1 + dot(&s1[0], &s2[0])
+             + dot(&s2[0], &s3[0])
              + dot(&s3[0], &s1[0])
              + I * dot(&s1[0], &crossp[0])
            ) / rho;
@@ -193,23 +227,84 @@ double compute_BergLuscher_angle(double *s1, double *s2, double *s3){
 double skyrmion_number_BergLuscher(double *spin, double *charge,
                        int nx, int ny, int nz, int *ngbs) {
 
-    int n = nx * ny * nz;
-    int i, spin_index;
-    double total_sum = 0;
+    /* Compute the topological charge (or skyrmion number) by adding triangles
+     * of neighbouring spins for every lattice site, which cover triangle areas
+     * in a unit sphere (i.e. we map the lattice area into a unit sphere
+     * surface, using a triangulation). The exponential that defines every
+     * spherical angle was firstly mentioned by Berg and Luscher [Nucl Phys B
+     * 190, 412 (1981)] for a discrete square lattice but we generalise it here
+     * for hexagonal crystals.
+     *
+     * NEIGHBOURS DEFINITION:
+     *
+     * *ngbs is the array with the neighbours information for every lattice
+     * site. For cuboid meshes, the array is like:
+     *
+     *  NN:      j =0  j=1 ...                 j=0  j=1  ...
+     *         [ 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, 1-x, 1+x, 1-y, ...  ]
+     *  spin:   i=6 * 0                        i=6 * 1
+     * ...
+     *
+     * where  0-y  is the index of the neighbour of the 0th spin, in the -y
+     * direction, for example, so the neighbours of the i-th spin start at the
+     * (6 * i) position of the array. This is similar for hexagonal meshes.
+     *
+     * For a cuboid mesh,  for the i-th spin, we generate the nearest ngbs
+     * triangles using the triangles: [i, j=1, j=3] , [i, j=0, j=2]
+     *                                     +x   +y         -x   -y
+     * i.e. we cover the top right and bottom left areas ( we could also use
+     * the top left and bottom right)
+     *
+     * For a hexagonal mesh,  for the i-th spin, we generate the nearest ngbs
+     * triangles using the triangles: [i, j=0, j=2] , [i, j=0, j=2]
+     *                                     E    NE         W    SW
+     *                                   (East) ...
+     * whose area covers a unit cell in a hexagonal crystal arrangement.
+     *
+     *
+     * Since the neighbours agree in indexes we can use the same function
+     * for both meshes.
+     *
+     * ------------------------------------------------------------------------
+     *
+     * Having the triangles defined, we compute the spherical triangle area
+     * spanned by the three spins using the compute_BergLuscher_angle function.
+     *
+     * The total topological charge Q is computed summing all triangles:
+     *                 __
+     *         Q  =   \     1   [ Omega(S_i, S_0, S_2) + Omega(S_i, S_1, S_3) ]
+     *                /__  ---
+     *                    4 PI
+     *                 i
+     *
+     * ehich are normalised by 4 PI, so we get an integer number for the number
+     * of times the unit sphere is covered by the spin directions in every
+     * triangle, e.g. if we have two skyrmions we get approximately Q = 2.
+     *
+     */
+
+    int n = nx * ny * nz; int i, spin_index; double total_sum = 0;
 
+    // Sweep through every lattice site
     for(i = 0; i < n; i++){
 
+        // x-y-z components for the i-th spin start at:
         spin_index = 3 * i;
 
+        // Reset the charge array
         charge[i] = 0;
 
+        // Compute the spherical triangle area for the triangle formed
+        // by the i-th spin and neighbours 0 and 2, i.e.
+        // [i j=0 j=2]. First check that the NNs exist:
         if(ngbs[6 * i] > 0 && ngbs[6 * i + 2] > 0){
             charge[i] += compute_BergLuscher_angle(&spin[spin_index],
                                                    &spin[3 * ngbs[6 * i]],
                                                    &spin[3 * ngbs[6 * i + 2]]
                                                    );
         }
 
+        // Triangle: [i j=1 j=3]
         if(ngbs[6 * i + 1] > 0 && ngbs[6 * i + 3] > 0){
             charge[i] += compute_BergLuscher_angle(&spin[spin_index],
                                                    &spin[3 * ngbs[6 * i + 1]],

fidimag/atomistic/llg.py

@@ -289,20 +289,41 @@ def compute_energy(self):
 
         return energy
 
-    def skyrmion_number(self):
-        nx = self.mesh.nx
-        ny = self.mesh.ny
-        nz = self.mesh.nz
-        number = clib.compute_skyrmion_number(
-            self.spin, self._skx_number, nx, ny, nz, self.mesh.neighbours)
-        return number
+    def skyrmion_number(self, method='FiniteSpinChirality'):
+        """
+        Calculate the skyrmion number using different methods:
+
+        method      :: 'FiniteSpinChirality', 'BergLuscher'
+
+        which are specifically defined for discrete spin lattices.
 
-    def skyrmion_number_BergLuscher(self):
+        Calling this function will fill the _skx_number array with the skyrmion
+        number density per lattice site.
+
+        See the corresponding C code at fidimag/atomistic/lib/util.c for a
+        detailed documentation about the methods.
+
+        """
         nx = self.mesh.nx
         ny = self.mesh.ny
         nz = self.mesh.nz
-        number = clib.compute_skyrmion_number_BergLuscher(
-            self.spin, self._skx_number, nx, ny, nz, self.mesh.neighbours)
+
+        if method == 'FiniteSpinChirality':
+            if self.mesh.mesh_type == 'cuboid':
+                number = clib.compute_skyrmion_number(
+                    self.spin, self._skx_number, nx, ny, nz,
+                    self.mesh.neighbours)
+            else:
+                raise ValueError('FiniteSpinChirality method only'
+                                 ' defined for cuboid meshes')
+
+        elif method == 'BergLuscher':
+            number = clib.compute_skyrmion_number_BergLuscher(
+                self.spin, self._skx_number, nx, ny, nz,
+                self.mesh.neighbours)
+        else:
+            raise ValueError('Specify a valid method')
+
         return number
 
     def spin_at(self, i, j, k):