add a rkky-style exchange for micro code

Author: Weiwei Wang

examples/micromagnetic/PRB_88_184422/fidimag/prb.py

@@ -32,7 +32,7 @@ def test_prb88_184422():
 
     sim.set_m((0, 0, 1))
 
-    sim.add(UniformExchange(A))
+    sim.add(UniformExchange(A=A))
     sim.add(DMI(-D, type='interfacial'))
     sim.add(UniaxialAnisotropy(K, axis=(0, 0, 1)))
 

examples/micromagnetic/dw_stt/main.py

@@ -3,12 +3,12 @@
 import matplotlib.pyplot as plt
 
 import numpy as np
-from micro import Sim
-from common import CuboidMesh
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
 
 # The energies, we can use DMI in a future simulation
-from micro import UniformExchange
-from micro import UniaxialAnisotropy
+from fidimag.micro import UniformExchange
+from fidimag.micro import UniaxialAnisotropy
 # from micro import DMI
 
 mu0 = 4 * np.pi * 1e-7
@@ -21,11 +21,11 @@ def init_m(pos):
     x = pos[0]
 
     if x < 400:
-        return (1, 0, 0)
+        return (-1, 0, 0)
     elif 400 <= x < 500:
         return (0, 1, 1)
     else:
-        return (-1, 0, 0)
+        return (1, 0, 0)
 
 
 def relax_system(mesh):
@@ -38,7 +38,7 @@ def relax_system(mesh):
     sim.driver.alpha = 0.5
     sim.driver.gamma = 2.211e5
     sim.Ms = 8.6e5
-    sim.do_precession = False
+    sim.driver.do_precession = False
 
     # The initial state passed as a function
     sim.set_m(init_m)
@@ -59,8 +59,8 @@ def relax_system(mesh):
     # ONE_DEGREE_PER_NS = 17453292.52
 
     # Start relaxation and save the state in m0.npy
-    sim.relax(dt=1e-14, stopping_dmdt=0.00001, max_steps=5000,
-              save_m_steps=None, save_vtk_steps=None)
+    sim.relax(dt=1e-14, stopping_dmdt=0.01, max_steps=5000,
+              save_m_steps=None, save_vtk_steps=50)
 
     np.save('m0.npy', sim.spin)
 
@@ -106,8 +106,7 @@ def deal_plot(_list, output_file):
 
 
 # THIS NEEDS REVISION
-def deal_plot_dynamics(spin, m_component='mz',
-                       output_file):
+def deal_plot_dynamics(spin, m_component='mz', output_file='output.pdf'):
     """
     The dynamics of the m_component of the i-th spin
     of the chain, from the magnetisation snapshot
@@ -140,10 +139,10 @@ def deal_plot_dynamics(spin, m_component='mz',
 def excite_system(mesh):
 
     # Specify the stt dynamics in the simulation
-    sim = Sim(mesh, name='dyn', driver='llg_stt')
+    sim = Sim(mesh, name='dyn2', driver='llg_stt')
 
     sim.driver.set_tols(rtol=1e-12, atol=1e-14)
-    sim.driver.alpha = 0.05
+    sim.driver.alpha = 0.2
     sim.driver.gamma = 2.211e5
     sim.Ms = 8.6e5
 
@@ -163,16 +162,16 @@ def excite_system(mesh):
 
     # Set the current in the x direction, in A / m
     # beta is the parameter in the STT torque
-    sim.jx = -1e12
-    sim.beta = 1
+    sim.driver.jx = -1e12
+    sim.driver.beta = 0.01
 
     # The simulation will run for 5 ns and save
     # 500 snapshots of the system in the process
     ts = np.linspace(0, 5e-9, 501)
 
     for t in ts:
         print 'time', t
-        sim.run_until(t)
+        sim.driver.run_until(t)
         sim.save_vtk()
         sim.save_m()
 

examples/micromagnetic/vortex/main.py

@@ -3,10 +3,10 @@
 import matplotlib.pyplot as plt
 
 import numpy as np
-from micro import Sim
-from common import CuboidMesh
-from micro import UniformExchange, Demag
-from micro import Zeeman, TimeZeeman
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
+from fidimag.micro import UniformExchange, Demag
+from fidimag.micro import Zeeman, TimeZeeman
 from fidimag.common.fileio import DataReader
 
 mu0 = 4 * np.pi * 1e-7
@@ -46,7 +46,7 @@ def relax_system(mesh):
     sim.driver.alpha = 0.5
     sim.driver.gamma = 2.211e5
     sim.Ms = spatial_Ms
-    sim.do_precession = False
+    sim.driver.do_precession = False
 
     sim.set_m(init_m)
     # sim.set_m(np.load('m0.npy'))
@@ -55,7 +55,7 @@ def relax_system(mesh):
     exch = UniformExchange(A=A)
     sim.add(exch)
 
-    demag = Demag()
+    demag = Demag(pbc_2d=True)
     sim.add(demag)
 
     mT = 795.7747154594767
@@ -84,7 +84,7 @@ def excite_system(mesh):
     exch = UniformExchange(A=A)
     sim.add(exch)
 
-    demag = Demag()
+    demag = Demag(pbc_2d=True)
     sim.add(demag)
 
     mT = 795.7747154594767
@@ -110,9 +110,9 @@ def gaussian_fun(t):
 
     mesh = CuboidMesh(nx=80, ny=80, nz=2, dx=2.5, dy=2.5, dz=5.0, unit_length=1e-9)
 
-    # relax_system(mesh)
+    relax_system(mesh)
 
-    excite_system(mesh)
+    #excite_system(mesh)
 
     # apply_field1(mesh)
     # deal_plot()

fidimag/micro/__init__.py

@@ -1,8 +1,10 @@
 from .sim import Sim
 from .demag import Demag
 from .exchange import UniformExchange
+from .exchange_rkky import ExchangeRKKY
 from .zeeman import Zeeman, TimeZeeman
 from .anisotropy import UniaxialAnisotropy
 from .dmi import DMI
 from .baryakhtar import LLBar, LLBarFull
 from .simple_demag import SimpleDemag
+

fidimag/micro/exchange_rkky.py

@@ -0,0 +1,61 @@
+import fidimag.extensions.micro_clib as micro_clib
+from .energy import Energy
+#from constant import mu_0
+
+
+class ExchangeRKKY(Energy):
+
+    """
+    RKKYExchange(sigma, Delta, z_bottom=0, z_top = -1, name='RKKYExchange')
+
+    Compute the RKKY-style exchange interaction defined by
+
+    E = sigma/Delta *(1-m_t * m_b)
+
+    where E is the energy density, sigma is the surface exchange coefficient between the two surfaces, 
+    Delata is the space thickness. m_t and m_b are the unit vectors of the top and bottom layers, respectively.
+    
+    Inputs:
+        sigma: float
+            sigma is the surface exchange stiffness constant.
+        Delta: float
+            Delta is the the space thickness in Meter.
+        z_bottom: int
+            z_bottom is the index of the bottom layer
+        z_top: int
+            z_top is the index of the top layer.
+    
+    """
+
+    def __init__(self, sigma, Delta=1e-9, z_bottom=0, z_top = -1, name='RKKYExchange'):
+        self.sigma = sigma/Delta
+        self.name = name
+        self.z_bottom = z_bottom
+        self.z_top = z_top
+        self.jac = True
+
+    def setup(self, mesh, spin, Ms):
+        super(ExchangeRKKY, self).setup(mesh, spin, Ms)
+        if self.z_top<0:
+            self.z_top+= self.nz
+ 
+
+    def compute_field(self, t=0, spin=None):
+        if spin is not None:
+            m = spin
+        else:
+            m = self.spin
+
+        micro_clib.compute_exchange_field_micro_rkky(m,
+                                                self.field,
+                                                self.energy,
+                                                self.Ms_inv,
+                                                self.sigma,
+                                                self.nx,
+                                                self.ny,
+                                                self.nz,
+                                                self.z_bottom,
+                                                self.z_top
+                                                )
+
+        return self.field

fidimag/micro/lib/exch.c

@@ -10,22 +10,22 @@ void compute_exch_field_micro(double *m, double *field, double *energy,
      *
      * Ms_inv     :: Array with the (1 / Ms) values for every mesh node.
      *               The values are zero for points with Ms = 0 (no material)
-     *  
+     *
      * A          :: Exchange constant
-     * 
+     *
      * dx, dy, dz :: Mesh spacings in the corresponding directions
-     * 
+     *
      * n          :: Number of mesh nodes
      *
      * ngbs       :: The array of neighbouring spins, which has (6 * n)
-     *               entries. Specifically, it contains the indexes of 
+     *               entries. Specifically, it contains the indexes of
      *               the neighbours of every mesh node, in the following order:
      *                      -x, +x, -y, +y, -z, +z
-     *  
+     *
      *               Thus, 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. The index value for a
      *              neighbour where Ms = 0, is evaluated as -1. The array
@@ -48,15 +48,15 @@ void compute_exch_field_micro(double *m, double *field, double *energy,
      *              since it is the left neighbour which is the PBC in x, etc..)
      *
      *  For the exchange computation, the field is defined as:
-     *          H_ex = (2 * A / (mu0 * Ms)) * nabla^2 (mx, my, mz) 
+     *          H_ex = (2 * A / (mu0 * Ms)) * nabla^2 (mx, my, mz)
      *
      *  Therefore, for the i-th mesh node (spin), we approximate the
      *  derivatives as:
-     *          nabla^2 mx = (1 / dx^2) * ( m[i-x] - 2 * m[i] + m[i+x] ) + 
-     *                       (1 / dy^2) * ( m[i-y] - 2 * m[i] + m[i+y] ) +  
+     *          nabla^2 mx = (1 / dx^2) * ( m[i-x] - 2 * m[i] + m[i+x] ) +
+     *                       (1 / dy^2) * ( m[i-y] - 2 * m[i] + m[i+y] ) +
      *                       (1 / dz^2) * ( m[i-z] - 2 * m[i] + m[i+z] )
-     * 
-     *  Where i-x is the neighbour in the -x direction. This is similar 
+     *
+     *  Where i-x is the neighbour in the -x direction. This is similar
      *  for my and mz.
      *  We can notice that the sum is the same if we do:
      *        ( m[i-x] - m[i] ) + ( m[i+x] - m[i] )
@@ -101,7 +101,7 @@ void compute_exch_field_micro(double *m, double *field, double *energy,
 	        field[3 * i + 2] = 0;
 	        continue;
 	    }
-        
+
         /* Here we iterate through the neighbours */
         for (int j = 0; j < 6; j++) {
             /* Remember that index=-1 is for sites without material */
@@ -113,15 +113,15 @@ void compute_exch_field_micro(double *m, double *field, double *energy,
                 /* Check that the magnetisation of the neighbouring spin
                  * is larger than zero */
                 if (Ms_inv[ngbs[idn + j]] > 0){
-                
-                    /* Neighbours in the -x and +x directions 
+
+                    /* Neighbours in the -x and +x directions
                      * giving: ( m[i-x] - m[i] ) + ( m[i+x] - m[i] )
                      * when ngbs[idn + j] > 0 for j = 0 and j=1
                      * If, for example, there is no
                      * neighbour at -x (j=0) in the 0th node (no PBCs),
                      * the second derivative would only be avaluated as:
                      *      (1 / dx * dx) * ( m[i+x] - m[i] )
-                     * which, according to 
+                     * which, according to
                      * [M.J. Donahue and D.G. Porter; Physica B, 343, 177-183 (2004)]
                      * when performing the integration of the energy, we still
                      * have error of the order O(dx^2)
@@ -160,3 +160,77 @@ void compute_exch_field_micro(double *m, double *field, double *energy,
         field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV;
     }
 }
+
+inline int get_index(int nx, int ny, int i, int j, int k){
+ return k * nx*ny + j * nx + i;
+}
+
+void compute_exch_field_rkky_micro(double *m, double *field, double *energy, double *Ms_inv,
+                         double sigma, int nx, double ny, double nz, int z_bottom, int z_top){
+
+    /* Compute the micromagnetic exchange field and energy using the
+     * matrix of neighbouring spins and a second order approximation
+     * for the derivative
+     *
+     * Ms_inv     :: Array with the (1 / Ms) values for every mesh node.
+     *               The values are zero for points with Ms = 0 (no material)
+     *
+     * sigma          :: Exchange constant
+     *
+     * nx, ny, nz :: Mesh dimensions.
+     *  The exchange field at the top (bottom) layer can be computed as:
+     *
+     *          H_top = (sigma / (mu0 * Ms)) * m_bottom
+     *          H_bottom = (sigma / (mu0 * Ms)) * m_top
+     *
+     *  The *m array contains the spins as:
+     *      [mx0, my0, mz0, mx1, my1, mz1, mx2, ...]
+     *  so if we want the starting position of the magnetisation for the
+     *  i-th spin, we only have to do (3 * i) for mx, (3 * i + 1) for my
+     *  and (3 * i + 2) for mz
+     *
+     *
+     *
+     */
+
+		 int n = nx*ny*nz;
+		 for (int i = 0; i < n; i++){
+        energy[i] = 0;
+				field[3*i]=0;
+				field[3*i+1]=0;
+				field[3*i+2]=0;
+		 }
+
+		 #pragma omp parallel for
+		 for (int i = 0; i < nx; i++) {
+			 for (int j = 0; j < ny; j++){
+				 double mtx=0, mty=0, mtz=0;
+				 double mbx=0, mby=0, mbz=0;
+				 int id1 = get_index(nx,ny, i, j, z_bottom);
+				 int id2 = get_index(nx,ny, i, j, z_top);
+				 mtx = m[3*id2];
+				 mty = m[3*id2+1];
+				 mtz = m[3*id2+2];
+
+				 mbx = m[3*id1];
+				 mby = m[3*id1+1];
+				 mbz = m[3*id1+2];
+
+				 if (Ms_inv[id1] != 0.0){
+					 energy[id1]  = sigma*(1-mtx*mbx-mty*mby-mtz*mbz);
+					 field[3*id1]   = sigma * mtx * Ms_inv[id1] * MU0_INV;
+					 field[3*id1+1] = sigma * mty * Ms_inv[id1] * MU0_INV;
+					 field[3*id1+2] = sigma * mtz * Ms_inv[id1] * MU0_INV;
+				 }
+
+				 if (Ms_inv[id2] != 0.0){
+					 energy[id2]  = sigma*(1-mtx*mbx-mty*mby-mtz*mbz);
+					 field[3*id2]   = sigma * mbx * Ms_inv[id2] * MU0_INV;
+					 field[3*id2+1] = sigma * mby * Ms_inv[id2] * MU0_INV;
+					 field[3*id2+2] = sigma * mbz * Ms_inv[id2] * MU0_INV;
+				 }
+			 }
+		 }
+
+
+}

fidimag/micro/lib/micro_clib.h

@@ -16,6 +16,9 @@ inline double cross_z(double a0, double a1, double a2,
 void compute_exch_field_micro(double *m, double *field, double *energy, double *Ms_inv,
                          double A, double dx, double dy, double dz, int n, int *ngbs);
 
+void compute_exch_field_rkky_micro(double *m, double *field, double *energy, double *Ms_inv,
+                         double sigma, int nx, double ny, double nz, int z_bottom, int z_top);
+
 void dmi_field(double *m, double *field, double *energy, double *Ms_inv,
                double *D, double dmi_vector[18], int n_dmi_ngbs,
                double dx, double dy, double dz,