Added micromagnetic calculation of the skyrmion number, discretising the linear derivatives of the magnetisation for the calculation of the topological charge. Changed the *type* argument for the DMI in the micromagnetic library, to *dmi_type* as in the atomistic library. This will break some docs or notebook examples. Added test for the micromagnetic sk number and documented the atomistic test

Author: davidcortesortuno

fidimag/micro/dmi.py

@@ -12,14 +12,14 @@ class DMI(Energy):
         compute the DMI field in micromagnetics
     """
 
-    def __init__(self, D, name='dmi', type='bulk'):
+    def __init__(self, D, name='dmi', dmi_type='bulk'):
         """
         type could be 'interfacial' or 'bulk'
         """
         self.D = D
         self.name = name
         self.jac = True
-        self.type = type
+        self.dmi_type = dmi_type
 
     def setup(self, mesh, spin, Ms):
         super(DMI, self).setup(mesh, spin, Ms)
@@ -32,7 +32,7 @@ def compute_field(self, t=0, spin=None):
         else:
             m = self.spin
 
-        if self.type == 'bulk':
+        if self.dmi_type == 'bulk':
             micro_clib.compute_dmi_field_bulk(m,
                                               self.field,
                                               self.energy,
@@ -45,7 +45,7 @@ def compute_field(self, t=0, spin=None):
                                               self.neighbours
                                               )
 
-        elif self.type == 'interfacial':
+        elif self.dmi_type == 'interfacial':
             micro_clib.compute_dmi_field_interfacial(m,
                                                      self.field,
                                                      self.energy,
@@ -59,6 +59,6 @@ def compute_field(self, t=0, spin=None):
                                                      )
         else:
             raise Exception(
-                "Unsppourted dmi type:{}, avaiable type: 'bulk','interfacial'.".format(self.type))
+                "Unsppourted dmi type:{}, avaiable type: 'bulk','interfacial'.".format(self.dmi_type))
 
         return self.field

fidimag/micro/lib/micro_clib.h

@@ -26,3 +26,6 @@ void dmi_field_interfacial(double *m, double *field, double *energy, double *Ms_
 
 void compute_uniaxial_anis(double *m, double *field, double *energy, double *Ms_inv, 
 	double *Ku, double *axis, int nx, int ny, int nz);
+
+double skyrmion_number(double *spin, double *charge,
+                       int nx, int ny, int nz, int *ngbs);

fidimag/micro/lib/micro_clib.pyx

@@ -22,7 +22,9 @@ cdef extern from "micro_clib.h":
                                double *energy, double *Ms_inv, 
                                double *Ku, double *axis,
                                int nx, int ny, int nz)
-
+    
+    double skyrmion_number(double *m, double *charge,
+                           int nx, int ny, int nz, int *ngbs)
 
 
 def compute_exchange_field_micro(np.ndarray[double, ndim=1, mode="c"] m,
@@ -70,3 +72,10 @@ def compute_anisotropy_micro(np.ndarray[double, ndim=1, mode="c"] m,
 
     compute_uniaxial_anis(&m[0], &field[0], &energy[0], &Ms_inv[0],
                           &Ku[0], &axis[0], nx, ny, nz)
+
+def compute_skyrmion_number(np.ndarray[double, ndim=1, mode="c"] m,
+                            np.ndarray[double, ndim=1, mode="c"] charge,
+                            nx, ny, nz,
+            		    np.ndarray[int, ndim=2, mode="c"] ngbs):
+
+    return skyrmion_number(&m[0], &charge[0], nx, ny, nz, &ngbs[0, 0])

fidimag/micro/lib/util.c

@@ -0,0 +1,101 @@
+#include "micro_clib.h"
+
+// Compute: S \cdot (S_i \times S_j)
+inline double volume(double S[3], double Si[3], double Sj[3]) {
+	double tx = S[0] * (-Si[2] * Sj[1] + Si[1] * Sj[2]);
+	double ty = S[1] * (Si[2] * Sj[0] - Si[0] * Sj[2]);
+	double tz = S[2] * (-Si[1] * Sj[0] + Si[0] * Sj[1]);
+	return tx + ty + tz;
+}
+
+double skyrmion_number(double *spin, double *charge,
+                       int nx, int ny, int nz, int *ngbs) {
+
+    /* Compute the skyrmion number Q, defined as:
+     *                      _
+     *           1         /       dM     dM
+     *   Q   =  ---  *    /   M .  --  X  --   dx dy
+     *          4 PI   _ /         dx     dy
+     *
+     *
+     * for a two dimensional layer (it can be a slice of a
+     * bulk system)
+     *
+     * Therefore, a finite differences discretisation of the
+     * continuum magnetisation field, using central
+     * differences, leads to:
+     *
+     * Q =   M_i \cdot ( M_{i+1} \times M_{j+1} )
+     *       + M_i \cdot ( M_{i-1} \times M_{j-1} )
+     *       - M_i \cdot ( M_{i-1} \times M_{j+1} )
+     *       - M_i \cdot ( M_{i+1} \times M_{j-1} )
+     *
+     *
+     * INPUTS:
+     *
+     * 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
+     */
+
+	int i, j;
+	int index, id;
+
+	double sum = 0;
+
+	double S[3];
+
+	int nxy = nx * ny;
+
+    /* Store the spin directions of the nearest neighbours
+     * in the order: [-x +x -y +y]
+     */ 
+    double S_nn[12];
+
+	for (i = 0; i < nxy; i++) {
+        index = 3 * 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];
+        S[2] = spin[index + 2];
+
+        for (j = 0; j < 4; j++) {
+            if (ngbs[id_nn + j] > 0) {
+                id = 3 * ngbs[id_nn + j];
+                S_nn[3 * j    ] = spin[id];
+                S_nn[3 * j + 1] = spin[id + 1];
+                S_nn[3 * j + 2] = spin[id + 2];
+            }
+        }
+
+        charge[i] =   volume(S, &S_nn[3], &S_nn[9])
+                    + volume(S, &S_nn[0], &S_nn[6])
+                    - volume(S, &S_nn[0], &S_nn[9])
+                    - volume(S, &S_nn[3], &S_nn[6]);
+
+        charge[i] /= (16 * WIDE_PI);
+
+        sum += charge[i];
+    }
+
+	return sum;
+
+}

fidimag/micro/llg.py

@@ -3,6 +3,7 @@
 import os
 import time
 import fidimag.extensions.clib as clib
+import fidimag.extensions.micro_clib as micro_clib
 import numpy as np
 from fidimag.common.fileio import DataSaver, DataReader
 from fidimag.common.save_vtk import SaveVTK
@@ -303,7 +304,7 @@ def skyrmion_number(self):
         nx = self.mesh.nx
         ny = self.mesh.ny
         nz = self.mesh.nz
-        number = clib.compute_skyrmion_number(
+        number = micro_clib.compute_skyrmion_number(
             self.spin, self._skx_number, nx, ny, nz, self.mesh.neighbours)
         return number
 

tests/test_oommf_without_run.py

@@ -141,7 +141,7 @@ def test_dmi_field_oommf(D=4.1e-3, Ms=2.6e5):
     sim = Sim(mesh)
     sim.Ms = Ms
 
-    dmi = DMI(D=D, type='interfacial')
+    dmi = DMI(D=D, dmi_type='interfacial')
     sim.add(dmi)
 
     def init_m(pos):

tests/test_prb88_184422.py

@@ -35,7 +35,7 @@ def run_fidimag(mesh):
     sim.set_m((0, 0, 1))
 
     sim.add(UniformExchange(A))
-    sim.add(DMI(D, type='interfacial'))
+    sim.add(DMI(D, dmi_type='interfacial'))
     sim.add(UniaxialAnisotropy(K, axis=(0, 0, 1)))
 
     sim.relax(dt=1e-13, stopping_dmdt=0.01, max_steps=5000,