add cubic anisotropy for atomistic model

Author: Weiwei Wang

doc/core_eqs.rst

@@ -110,7 +110,21 @@ corresponding field is
 .. math::
    \vec{H}_{i,\text{an}} = \frac{2 \mathcal{K}_{u}}{\mu_s} \left(\vec{S}_{i}\cdot\hat{u}\right)\hat{u}
 
-|
+The Hamiltonian for the cubic anisotropy is given by
+
+.. math::
+   \mathcal{H}_{\text{an}}^c = \mathcal{K}_{c} \sum_i (S_x^4+S_y^4+S_z^4)
+
+which is equivalent to the popular form
+
+.. math::
+   \mathcal{H}_{\text{an}}^c = -2 \mathcal{K}_{c} \sum_i (S_x^2 S_y^2 + S_y^2 S_z^2 + S_z^2 S_x^2)
+
+The effective fields thus can be computed as
+
+.. math::
+   \vec{H}_{i,\text{an}}^c = -\frac{4 \mathcal{K}_{c}}{\mu_s} \left ( S_x^3 \hat{x} + S_y^3 \hat{y} + S_z^3 \hat{z} \right)
+
 
 In micromagnetics, the uniaxial anisotropy energy of the system is defined as
 

fidimag/atomistic/__init__.py

@@ -1,6 +1,6 @@
 from .sim import Sim
 from .exchange import UniformExchange, Exchange
-from .anisotropy import Anisotropy
+from .anisotropy import Anisotropy, CubicAnisotropy
 from .zeeman import Zeeman, TimeZeeman
 from .demag import Demag
 from .demag_hexagonal import DemagHexagonal

fidimag/atomistic/anisotropy.py

@@ -72,3 +72,33 @@ def compute_field(self, t=0, spin=None):
                                 )
 
         return self.field * self.mu_s_inv
+
+
+class CubicAnisotropy(Energy):
+    """
+    Compute the Cubic Anisotropy, see documentation for detailed equations.
+    """
+
+    def __init__(self, Kc, name='CubicAnisotropy'):
+        self.Kc = Kc
+        self.name = name
+        self.jac = True
+
+
+    def setup(self, mesh, spin, mu_s):
+        super(CubicAnisotropy, self).setup(mesh, spin, mu_s)
+        self._Kc = helper.init_scalar(self.Kc, self.mesh)
+
+    def compute_field(self, t=0, spin=None):
+        if spin is not None:
+            m = spin
+        else:
+            m = self.spin
+
+        clib.compute_anisotropy_cubic(m,
+                                self.field,
+                                self.energy,
+                                self._Kc,
+                                self.n)
+
+        return self.field * self.mu_s_inv

fidimag/atomistic/lib/anis.c

@@ -3,8 +3,8 @@
 
 void compute_anis(double *spin, double *field, double *energy,
 	double *Ku, double *axis, int n) {
-	
-    /* Remember that the magnetisation order is 
+
+    /* Remember that the magnetisation order is
      *      mx1, my1, mz1, mx2, my2, mz2, mx3,...
      * so we get the corresponding components multiplying
      * by 3 in every iteration. The anisotropy *axis* has the
@@ -20,8 +20,8 @@ void compute_anis(double *spin, double *field, double *energy,
         double m_u = (spin[3 * i] * axis[3 * i] +
                       spin[3 * i + 1] * axis[3 * i + 1] +
                       spin[3 * i + 2] * axis[3 * i + 2]);
-    
-    
+
+
 		field[3 * i]     = 2 * Ku[i] * m_u * axis[3 * i]    ;
 		field[3 * i + 1] = 2 * Ku[i] * m_u * axis[3 * i + 1];
 		field[3 * i + 2] = 2 * Ku[i] * m_u * axis[3 * i + 2];
@@ -31,3 +31,26 @@ void compute_anis(double *spin, double *field, double *energy,
 	}
 
 }
+
+
+void compute_anis_cubic(double *spin, double *field, double *energy,
+	double *Kc, int n) {
+
+    /* Remember that the magnetisation order is
+     *      mx1, my1, mz1, mx2, my2, mz2, mx3,...
+     * so we get the corresponding components multiplying
+     * by 3 in every iteration.
+     *
+     */
+      #pragma omp parallel for
+	    for (int i = 0; i < n; i++) {
+				int j = 3*i;
+				field[j] = - 4*Kc[i]*spin[j]*spin[j]*spin[j];
+				field[j+1] = - 4*Kc[i]*spin[j+1]*spin[j+1]*spin[j+1];
+				field[j+2] = - 4*Kc[i]*spin[j+2]*spin[j+2]*spin[j+2];
+
+	      energy[i] = -0.25*(field[j]*spin[j]+field[j+1]*spin[j+1]+field[j+2]*spin[j+2]) ;
+
+	    }
+
+}

fidimag/atomistic/lib/clib.h

@@ -36,6 +36,8 @@ double compute_exch_energy(double *spin, double Jx, double Jy, double Jz,
 
 void compute_anis(double *spin, double *field, double *energy, double *Ku,
                   double *axis, int n);
+void compute_anis_cubic(double *spin, double *field, double *energy,
+	          double *Kc, int n); 
 
 void dmi_field_bulk(double *spin, double *field, double *energy, double *D,
                     int *ngbs, int n);

fidimag/atomistic/lib/clib.pyx

@@ -64,6 +64,8 @@ cdef extern from "clib.h":
     void compute_anis(double *spin, double *field, double *energy,
                       double *Ku, double *axis, int n)
 
+    void compute_anis_cubic(double *spin, double *field, double *energy, double *Kc, int n)
+
     void llg_rhs(double * dm_dt, double * spin,
                  double *h, double *alpha, int *pins,
                  double gamma, int n, int do_precession, double default_c)
@@ -159,14 +161,11 @@ def compute_exchange_energy(np.ndarray[double, ndim=1, mode="c"] spin,
                                xperiodic, yperiodic)
 
 
-def compute_anisotropy(np.ndarray[double, ndim=1, mode="c"] spin,
-                        np.ndarray[double, ndim=1, mode="c"] field,
-                        np.ndarray[double, ndim=1, mode="c"] energy,
-                        np.ndarray[double, ndim=1, mode="c"] Ku,
-                        np.ndarray[double, ndim=1, mode="c"] axis,
-                        n):
-    compute_anis(&spin[0], &field[0], &energy[0], &Ku[0],
-                 &axis[0], n)
+def compute_anisotropy(double [:] spin, double [:] field, double [:] energy, double [:] Ku, double [:] axis, n):
+    compute_anis(&spin[0], &field[0], &energy[0], &Ku[0], &axis[0], n)
+
+def compute_anisotropy_cubic(double [:] spin, double [:] field, double [:] energy, double [:] Kc, n):
+    compute_anis_cubic(&spin[0], &field[0], &energy[0], &Kc[0], n)
 
 
 def compute_dmi_field(np.ndarray[double, ndim=1, mode="c"] spin,

fidimag/common/vtk.py

@@ -21,6 +21,7 @@ def __init__(self, mesh, header="", directory=".", filename="unnamed"):
             # for keyword argument dimensions: if the mesh is made up of
             # nx * ny * nz cells, it has (nx + 1) * (ny + 1) * (nz + 1)
             # vertices.
+            print(mesh.grid)
             structure = pyvtk.RectilinearGrid(* mesh.grid)
         else:
             raise NotImplementedError(

tests/test_anis.py

@@ -1,4 +1,4 @@
-from fidimag.atomistic import Anisotropy
+from fidimag.atomistic import Anisotropy, CubicAnisotropy
 from fidimag.common import CuboidMesh
 from fidimag.atomistic import Sim
 import numpy as np
@@ -17,6 +17,19 @@ def test_anis():
     field = anis.compute_field()
     assert field[0] == 2 * 99
 
+def test_anis_cubic():
+    mesh = CuboidMesh(nx=1, ny=1, nz=1)
+    spin = np.array([0.6,0.8,0])
+    anis = CubicAnisotropy(Kc=1.23)
+    mu_s = np.ones(3)
+    anis.setup(mesh, spin, mu_s)
+    field = anis.compute_field()
+    print field
+    assert  np.max(field-np.array([-1.06272,-2.51904, -0.]))<1e-6
+    energy = anis.compute_energy()
+    assert abs(energy-0.663216)<1e-5
+
 
 if __name__ == '__main__':
     test_anis()
+    test_anis_cubic()