add next-neighbour interaction for cuboid mesh in monte carlo method

Author: Weiwei Wang

fidimag/atomistic/lib/clib.h

@@ -90,7 +90,8 @@ void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha,
                   double *mu_s_inv, int *pins, double *eta, int n, double gamma, double dt);
 
 //======================================================================
-void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, double J, double D, double *h, double Kc, int n, double T, int hexagnoal_mesh);
+
+void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, double J, double J1, double D, double D1, double *h, double Kc, int n, double T, int hexagnoal_mesh);
 
 
 #endif

fidimag/atomistic/lib/clib.pyx

@@ -19,8 +19,8 @@ cdef extern from "time.h":
     time_t time(time_t *timer)
 
 cdef extern from "clib.h":
-    void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, double J, double D, double *h, double Kc, int n, double T, int hexagnoal_mesh)
-    
+    void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, 
+                    double J, double J1, double D, double D1, double *h, double Kc, int n, double T, int hexagnoal_mesh)
     double skyrmion_number(double *spin, double *charge,
                            int nx, int ny, int nz, int *ngbs)
 
@@ -340,8 +340,9 @@ cdef class monte_carlo(object):
     def run_step(self,np.ndarray[double, ndim=1, mode="c"] spin,
                 np.ndarray[double, ndim=1, mode="c"] new_spin,
                 np.ndarray[int, ndim=2, mode="c"] ngbs,
-                J, D, np.ndarray[double, ndim=1, mode="c"] h, 
+                np.ndarray[int, ndim=2, mode="c"] nngbs,
+                J, J1, D, D1, np.ndarray[double, ndim=1, mode="c"] h, 
                 Kc, n, T, hexagnoal_mesh):
 
-        run_step_mc(self._c_state, &spin[0], &new_spin[0], &ngbs[0,0], J, D, &h[0], Kc, n, T, hexagnoal_mesh)
+        run_step_mc(self._c_state, &spin[0], &new_spin[0], &ngbs[0,0], &nngbs[0,0], J, J1, D, D1, &h[0], Kc, n, T, hexagnoal_mesh)
 

fidimag/atomistic/lib/mc.c

@@ -97,29 +97,37 @@ double compute_deltaE_exchange_DMI_hexagnoal(double *spin, double *new_spin, int
 
 }
 
-double compute_deltaE_exchange_DMI(double *spin, double *new_spin, int *ngbs, double J, double D, int i){
+double compute_deltaE_exchange_DMI(double *spin, double *new_spin, int *ngbs, int *nngbs, double J, double J1, double D, double D1, int i){
 
     int id_nn = 6 * i;
     double energy1=0, energy2=0;
 
     for (int j = 0; j < 6; j++) {
         int k = ngbs[id_nn + j];
         if (k >= 0) {
-            energy1 -= J*dot(&spin[3*i], &spin[3*k]);//exchange energy
+            energy1 -= J*dot(&spin[3*i], &spin[3*k]) ;//exchange energy
             energy1 += D*dmi_energy_site(&spin[3*i], &spin[3*k], j); //DMI energy
 
             energy2 -= J*dot(&new_spin[3*i], &spin[3*k]);
             energy2 += D*dmi_energy_site(&new_spin[3*i], &spin[3*k], j);
+        }
+        k = nngbs[id_nn + j];
+        if (k >= 0) {
+            energy1 -= J1*dot(&spin[3*i], &spin[3*k]) ;//exchange energy
+            energy1 += D1*dmi_energy_site(&spin[3*i], &spin[3*k], j); //DMI energy
 
+            energy2 -= J1*dot(&new_spin[3*i], &spin[3*k]);
+            energy2 += D1*dmi_energy_site(&new_spin[3*i], &spin[3*k], j);
         }
+        
     }
 
     return energy2-energy1;
 
 }
 
 
-void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, double J, double D, double *h, double Kc, int n, double T, int hexagnoal_mesh){
+void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, double J, double J1, double D, double D1, double *h, double Kc, int n, double T, int hexagnoal_mesh){
 
     double delta_E, r;
     int update=0;
@@ -136,7 +144,7 @@ void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs
         if(hexagnoal_mesh){
            delta_E += compute_deltaE_exchange_DMI_hexagnoal(&spin[0], &new_spin[0], &ngbs[0], J, D, index);
         }else{
-           delta_E += compute_deltaE_exchange_DMI(&spin[0], &new_spin[0], &ngbs[0], J, D, index);
+           delta_E += compute_deltaE_exchange_DMI(&spin[0], &new_spin[0], &ngbs[0], &nngbs[0], J, J1, D, D1, index);
         }
 
 

fidimag/atomistic/monte_carlo.py

@@ -34,21 +34,27 @@ def __init__(self, mesh, name='unnamed'):
         self.hexagnoal_mesh = False
         if mesh.mesh_type == 'hexagonal':
             self.hexagnoal_mesh =  True
+            #FIX ME !!!!
+            self.nngbs = np.copy(mesh.neighbours)
+        else:
+            self.nngbs = mesh.next_neighbours
 
         self.step = 0
         self.skx_num = 0
         self.mc = clib.monte_carlo()
         self.set_options()
 
-    def set_options(self, J=50.0, D=0, Kc=0, H=None, seed=100, T=10.0, S=1):
+    def set_options(self, J=50.0, J1=0, D=0, D1=0, Kc=0, H=None, seed=100, T=10.0, S=1):
         """
         J, D and Kc in units of k_B
         H in units of Tesla.
         S is the spin length
         """
         self.mc.set_seed(seed)
         self.J = J
+        self.J1 = J1
         self.D = D
+        self.D1 = D1
         self.T = T
         self.Kc = Kc
         self.mu_s =  1.0
@@ -145,8 +151,8 @@ def run(self, steps=1000, save_m_steps=100, save_vtk_steps=100, save_data_steps=
 
         for step in range(1, steps + 1):
             self.step = step
-            self.mc.run_step(self.spin, self.random_spin, self.ngbs,
-                self.J, self.D, self._H, self.Kc, self.n, self.T, self.hexagnoal_mesh)
+            self.mc.run_step(self.spin, self.random_spin, self.ngbs, self.nngbs, 
+                self.J, self.J1, self.D, self.D1, self._H, self.Kc, self.n, self.T, self.hexagnoal_mesh)
             if save_data_steps is not None:
                 if step % save_data_steps == 0:
                     self.saver.save()