Merge pull request #45 from computationalmodelling/monte_carlo

Monte carlo

Author: Weiwei Wang

doc/index.rst

@@ -19,5 +19,6 @@ Contents:
    install
    core_eqs
    extended_eqs
+   monte_carlo
    develop
    ipynb/tutorial-basics

doc/monte_carlo.rst

@@ -0,0 +1,42 @@
+Monte Carlo Simulation
+=======================
+
+In the atomistic part of Fidimag, Monte Carlo based on Metropolis algorithm is integrated. 
+The supported interactions include the exchange interaction, the bulk DMI, the external field and
+cubic anisotropy. The total Hamiltonian of the system is therefore given by
+
+.. math::
+   \mathcal{H} = \mathcal{H}_{ex} + \mathcal{H}_{dmi} + \mathcal{H}_{ext} + \mathcal{H}_{c}
+
+where :math:`\mathcal{H}_{ex}` is the nearest-neighbor exchange interaction, 
+
+.. math::
+   \mathcal{H}_{ex} = -J \sum_{<i,j>}\vec{m}_i \cdot \vec{m}_j
+
+Note that the summation is taken only once for each pair and :math:`\vec{m}_i` is the unit vector of spin :math:`\vec{S}` at site :math:`i`. 
+  
+The Hamiltonian of bulk DMI can be expressed as, 
+
+.. math::
+   \mathcal{H}_{dmi}= \sum_{<i,j>} \vec{D}_{ij}\cdot [\vec{m}_i \times \vec{m}_j]
+
+where :math:`\vec{D}_{ij} = D \vec{e}_{ij}` with :math:`\vec{e}_{ij}` is the unit vector between :math:`\vec{S}_{i}` and :math:`\vec{S}_{j}`.
+
+The Hamiltonian of external field is
+
+.. math::
+   \mathcal{H}_{dmi}= - \sum_{i} \mu_s \vec{H}\cdot  \vec{m}_i
+
+where :math:`\mu_s = g \mu_B S`.
+
+The cubic anisotropy implemented in Fidimag is, 
+
+.. math::
+   \mathcal{H}_{c}= - \sum_{i} (K_c/2)  (m_{i,x}^4 + m_{i,y}^4 + m_{i,z}^4)
+
+which is equivalent to the form 
+
+.. math::
+   \mathcal{H}_{c}= \sum_{i} K_c (m_{i,x}^2 m_{i,y}^2 + m_{i,y}^2 m_{i,z}^2 + m_{i,z}^2 m_{i,x}^2)
+
+

examples/atomic/PRL_111_067203/single.py

@@ -105,4 +105,4 @@ def excite_system(mesh):
 if __name__ == '__main__':
     mesh = CuboidMesh(nx=150, ny=50, nz=1,  periodicity=(True, True, False))
     relax_system(mesh)
-    #excite_system(mesh)
+    excite_system(mesh)

examples/mc/skx.py

@@ -3,7 +3,7 @@
 import matplotlib.pyplot as plt
 
 import numpy as np
-from fidimag.atomistic.mc import MonteCarlo
+from fidimag.atomistic import MonteCarlo
 from fidimag.common import DataReader, CuboidMesh
 
 
@@ -19,9 +19,9 @@ def random_m(pos):
 def run(mesh):
 
     mc = MonteCarlo(mesh, name='test1')
-    mc.set_m(init_m)
-    mc.set_options(H=[0,0,0.0], J=50.0, D=0.27*50, T=5.0)
-    mc.run(steps=200000, save_m_steps=None, save_vtk_steps=1000, save_data_steps=10)
+    mc.set_m(random_m)
+    mc.set_options(H=[0,0,0.0], J=50.0, D=0.27*50, T=5.0, Kc=50*0.1)
+    mc.run(steps=20000, save_m_steps=None, save_vtk_steps=1000, save_data_steps=10)
 
 if __name__=='__main__':
 

fidimag/atomistic/__init__.py

@@ -4,7 +4,9 @@
 from .zeeman import Zeeman, TimeZeeman
 from .demag import Demag
 from .demag_hexagonal import DemagHexagonal
+from .hexagonal_mesh import HexagonalMesh
 from .demag_full import DemagFull
 from .dmi import DMI
+from .monte_carlo import MonteCarlo
 import fidimag.common.constant as const
 from .materials import UnitMaterial, Nickel

fidimag/atomistic/lib/clib.h

@@ -6,11 +6,9 @@
 #include<fftw3.h>
 //#include<omp.h>
 
-#define WIDE_PI 3.1415926535897932384626433832795L
+#include"fidimag_random.h"
 
-double single_random(void);
-void gauss_random_vec(double *x, int n);
-void random_spin_uniform(double *spin, int n);
+#define WIDE_PI 3.1415926535897932384626433832795L
 
 /* 3 components for the cross product calculations */
 inline double cross_x(double a0, double a1, double a2,
@@ -92,8 +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(double *spin, double *new_spin, int *ngbs, double J, double D, double *h, int n, double T);
 
+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

@@ -2,13 +2,25 @@ import numpy
 cimport numpy as np
 np.import_array()
 
-cdef extern from "clib.h":
-    void initial_random(int seed)
-    double single_random()
-    void gauss_random_vec(double *x, int n)
-    void random_spin_uniform(double *spin, int n)
-    void run_step_mc(double *spin, double *new_spin, int *ngbs, double J, double D, double *h, int n, double T)
+cdef extern from "fidimag_random.h":
+    ctypedef struct mt19937_state:
+        pass
+    
+    mt19937_state *create_mt19937_state()
+    void finalize_mt19937_state(mt19937_state *state)
+
+    void initial_rng_mt19973(mt19937_state *state, int seed)
+    double random_double_half_open(mt19937_state *state)
+    void gauss_random_vector(mt19937_state *state, double *x, int n)
+    void uniform_random_sphere(mt19937_state *state, double *spin, int n)
 
+cdef extern from "time.h":
+    ctypedef int time_t
+    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, 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)
 
@@ -264,31 +276,6 @@ def normalise_spin(np.ndarray[double, ndim=1, mode="c"] spin,
     normalise(&spin[0], &pins[0], n)
 
 
-def init_random(seed):
-    initial_random(seed)
-
-def random_number_array(np.ndarray[double, ndim=1, mode="c"] v):
-    
-    cdef int n = len(v)
-
-    gauss_random_vec(&v[0], n)
-
-def random_spin_uniform_sphere(np.ndarray[double, ndim=1, mode="c"] v, n):
-
-    random_spin_uniform(&v[0], n)
-
-def random_number():
-    cdef double res = single_random()
-    return res
-
-def run_mc_step(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,
-                n, T):
-    run_step_mc(&spin[0], &new_spin[0], &ngbs[0,0], J, D, &h[0], n, T)
-
-
 def compute_llg_rhs_dw(np.ndarray[double, ndim=1, mode="c"] dm,
                 np.ndarray[double, ndim=1, mode="c"] spin,
                 np.ndarray[double, ndim=1, mode="c"] field,
@@ -298,3 +285,64 @@ def compute_llg_rhs_dw(np.ndarray[double, ndim=1, mode="c"] dm,
                 np.ndarray[double, ndim=1, mode="c"] eta,
                 np.ndarray[int, ndim=1, mode="c"] pin, n, gamma, dt):
     llg_rhs_dw_c(&spin[0], &field[0], &dm[0], &T[0], &alpha[0],  &mu_s_inv[0], &pin[0], &eta[0], n, gamma, dt)
+
+
+cdef class rng_mt19937(object):
+    cdef mt19937_state *_c_state
+    cdef public int seed
+    def __init__(self, seed=None):
+        if seed:
+            self.seed = int(seed)
+        else:
+            self.seed = time(NULL)
+                
+        self._c_state = create_mt19937_state()
+        if self._c_state is NULL:
+            raise MemoryError()
+        
+        initial_rng_mt19973(self._c_state, self.seed)
+
+    def set_seed(self, seed):
+        self.seed = int(seed)
+        initial_rng_mt19973(self._c_state, self.seed)
+
+    def random(self):
+        """
+            return a random number in [0,1)
+        """
+        return random_double_half_open(self._c_state)
+    
+    def fill_vector_gaussian(self, np.ndarray[np.float64_t, ndim=1] vector):
+        gauss_random_vector(self._c_state, &vector[0], vector.shape[0])
+
+    def fill_vector_uniform_sphere(self, np.ndarray[double, ndim=1, mode="c"] spin, n):
+        uniform_random_sphere(self._c_state,&spin[0], n)
+
+    def __dealloc__(self):
+        if self._c_state is not NULL:
+            finalize_mt19937_state(self._c_state)
+            self._c_state = NULL
+
+cdef class monte_carlo(object):
+    cdef mt19937_state *_c_state            
+
+    def __init__(self, seed=-1):
+                
+        self._c_state = create_mt19937_state()
+        if self._c_state is NULL:
+            raise MemoryError()
+        
+        initial_rng_mt19973(self._c_state, seed)
+
+    def set_seed(self, seed):
+        initial_rng_mt19973(self._c_state, seed)
+ 
+    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,
+                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], &nngbs[0,0], J, J1, D, D1, &h[0], Kc, n, T, hexagnoal_mesh)
+

fidimag/atomistic/lib/fidimag_random.c

@@ -0,0 +1,145 @@
+#include <stdlib.h>  // rand(), srand()
+#include <time.h>
+#include <math.h>
+#include "fidimag_random.h"
+
+mt19937_state *create_mt19937_state(void) {
+    
+    mt19937_state *state = (mt19937_state*) malloc(sizeof(mt19937_state));
+    state->seed = 1;
+    state->index_t = 0;
+    state->matrix[0] = 0;
+    state->matrix[1] = 0x9908b0dfU;
+    
+    return state;
+}
+
+void finalize_mt19937_state(mt19937_state *state) {
+    free(state);
+}
+
+
+void initial_rng_mt19973(mt19937_state *state, int seed) {
+    
+    if (seed<0){
+        state->seed = (unsigned int) time(NULL);
+    }else{
+        state->seed = seed;
+    }
+    
+    state->MT[0] = state->seed & 0xFFFFFFFFU;
+    for (int i = 1; i < 624; i++) {
+        state->MT[i] = (state->MT[i - 1] ^ (state->MT[i - 1] >> 30)) + i;
+        state->MT[i] *= 0x6c078965U;
+        state->MT[i] &= 0xFFFFFFFFU;
+    }
+}
+
+//return a integer in [0, MT19973_RAND_MAX]
+inline unsigned int rand_int(mt19937_state *state) {
+    
+    unsigned int x;
+    int index_t = state->index_t;
+    
+    x = (state->MT[index_t] & 0x1U) + (state->MT[(index_t + 1) % 624] & 0xFFFFFFFEU);
+    state->MT[index_t] = (state->MT[(index_t + 397) % 624] ^ (x >> 1))^ state->matrix[x & 1];
+    
+    x = state->MT[index_t];
+    x ^= (x >> 11);
+    x ^= (x << 7) & 0x9d2c5680U;
+    x ^= (x << 15) & 0xefc60000U;
+    x ^= (x >> 18);
+    
+    state->index_t = (index_t + 1) % 624;
+    
+    return x;
+}
+
+//return a double number in (0,1) with uniform distribution
+double random_double_open(mt19937_state *state) {
+	return (((double) rand_int(state))+0.5) / 4294967296.0;
+}
+
+//return a double number in [0,1) with uniform distribution
+double random_double_half_open(mt19937_state *state) {
+    return ((double) rand_int(state)) / 4294967296.0;
+}
+
+//return a double number in [0,1] with uniform distribution
+inline double random_double(mt19937_state *state) {
+    return ((double) rand_int(state)) / 4294967295.0;
+}
+
+//return a integer in [0, n-1]
+int rand_int_n(mt19937_state *state, int n){
+    double x = (double) rand_int(state) / 4294967296.0;
+    return (int)(n*x);
+}
+
+const double inv_a[] = { -3.969683028665376e+01, 2.209460984245205e+02,
+		-2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01,
+		2.506628277459239e+00 };
+
+const double inv_b[] = { -5.447609879822406e+01, 1.615858368580409e+02,
+		-1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01 };
+
+const double inv_c[] = { -7.784894002430293e-03, -3.223964580411365e-01,
+		-2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00,
+		2.938163982698783e+00 };
+
+const double inv_d[] = { 7.784695709041462e-03, 3.224671290700398e-01,
+		2.445134137142996e+00, 3.754408661907416e+00 };
+
+//inverse normal distribution function, see http://home.online.no/~pjacklam/notes/invnorm/
+//it seems the link is broken,  see http://www.psychometrica.de/norming/NormalDistribution.java for the jave implementation
+// or see https://github.com/stan-dev/stan/issues/1157 for another example
+inline double invnorm(mt19937_state *state) {
+	double q, r;
+	double p = (((double) rand_int(state))+0.5) / 4294967296.0; //0<p<1
+
+    if (p < 0.02425) {
+		q = sqrt(-2 * log(p));
+		return (((((inv_c[0] * q + inv_c[1]) * q + inv_c[2]) * q + inv_c[3]) * q + inv_c[4]) * q
+				+ inv_c[5]) / ((((inv_d[0] * q + inv_d[1]) * q + inv_d[2]) * q + inv_d[3]) * q + 1);
+	} else if (p > 0.97575) {
+		q = sqrt(-2 * log(1 - p));
+		return -(((((inv_c[0] * q + inv_c[1]) * q + inv_c[2]) * q + inv_c[3]) * q + inv_c[4]) * q
+				+ inv_c[5]) / ((((inv_d[0] * q + inv_d[1]) * q + inv_d[2]) * q + inv_d[3]) * q + 1);
+	} else {
+		q = p - 0.5;
+		r = q * q;
+		return (((((inv_a[0] * r + inv_a[1]) * r + inv_a[2]) * r + inv_a[3]) * r + inv_a[4]) * r
+				+ inv_a[5]) * q / (((((inv_b[0] * r + inv_b[1]) * r + inv_b[2]) * r + inv_b[3]) * r
+				+ inv_b[4]) * r + 1);
+	}
+}
+
+void gauss_random_vector(mt19937_state *state, double *x, int n) {
+    for (int i = 0; i < n; i++) {
+        x[i] = invnorm(state);
+    }
+}
+
+//generate an array with numbers in [0, n-1]
+void random_integer_vector(mt19937_state *state, int *ids, int n) {
+    for (int i = 0; i < n; i++) {
+        ids[i] = rand_int_n(state, n);
+    }
+}
+
+
+/*
+ * generate an uniform distribution in a spherical surface
+ * n is the total spin number, so len(spin) == 3*n
+ */
+void uniform_random_sphere(mt19937_state *state, double *spin, int n){
+    for (int i=0;i<n;i++){
+        int j=3*i;
+        double phi= random_double(state)*2*WIDE_PI;
+        double ct = 2*random_double(state)-1;
+        double st = sqrt(1-ct*ct);
+        spin[j] = st*cos(phi); //mx = sin(theta)*cos(phi)
+        spin[j+1] = st*sin(phi); //my = sin(theta)*sin(phi)
+        spin[j+2] = ct; //mz = cos(theta)
+    }
+}