Completed MnSi example

docs/MnSi.rst

@@ -1,7 +1,9 @@
 MnSi
 ====
 
-This example will guide you through the experimental investigation conducted by Amato et al. and Dalmas de Réotier et. al reported in the following journal articles [Amato2014]_, [DalmasdeReotier2016]_.
+This example will guide you through the experimental investigation conducted 
+by Amato et al. and Dalmas de Réotier et. al reported in the following journal
+articles [Amato2014]_, [DalmasdeReotier2016]_.
 
 .. note::  This is an advanced example. It is assumed that you already
            familiarized with Muesr by following one of the other examples
@@ -32,18 +34,18 @@ visually checking the dipolar tensor elements for all the equivalent sites. This
 analytically but a numerical check can come in handy.
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 36-78
+   :lines: 37-79
    :emphasize-lines: 4, 7,11-18,21-24,27
-   :lineno-start: 35
+   :lineno-start: 37
    :language: python
 
-In line 38 we add one of the four symmetry equivalent positions for a 
+In line 40 we add one of the four symmetry equivalent positions for a 
 muon in the `4a` Wyckoff site and we let the code find the other sites
-in line 41 with the function `muon_find_equiv`.
-In line 45 we define arbitrary small Fourier components which are used
+in line 42 with the function `muon_find_equiv`.
+In line 46 we define arbitrary small Fourier components which are used
 to specify which atoms are magnetic (Mn in this case).
 We finally add this sort of magnetic order to the sample description in 
-lines 55-58. The value of the propagation vector given in line 57 can 
+lines 56-59 The value of the propagation vector given in line 58 can 
 be omitted (resulting in 0 as default )as it is not used anywhere 
 in this initial estimation of the dipolar tensor. 
 The non-zero values of the dipolar tensor, shown at standard output, are
@@ -79,9 +81,9 @@ In order to compare with the experiment, we will evaluate the dipolar
 tensor for 100 values along the 111 direction. 
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 83-97
+   :lines: 84-98
    :emphasize-lines: 3,6
-   :lineno-start: 83
+   :lineno-start: 84
    :language: python
 
 
@@ -102,23 +104,23 @@ In order to calculate the local fields at the muon site in the helical state, we
 Let us first create a few useful variables for the definition of the Fourier components (FC).
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 115-146
-   :lineno-start: 115
+   :lines: 116-147
+   :lineno-start: 116
    :language: python
 
 We can now define both a left handed and a right handed helix.
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 150-169
-   :lineno-start: 150
+   :lines: 151-171
+   :lineno-start: 151
    :language: python
 
 We finally add the muon positions and the two magnetic orders with the 
 commands
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 175-188
-   :lineno-start: 175
+   :lines: 177-190
+   :lineno-start: 177
    :language: python
 
 .. note:: When a new magnetic order is added, the index is automatically
@@ -131,20 +133,20 @@ muon site we use the function locfield function and specify the 'i' sum
 type.
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 190-195
-   :lineno-start: 190
+   :lines: 192-197
+   :lineno-start: 192
    :language: python
 
 By default, the contact coupling is 0 for all sites. In order to have a 
 contact term different from 0 we have to set the parameter ACont
 for all the muon sites that we defined.
 
 .. literalinclude:: ../examples/MnSi/run_example.py
-   :lines: 198-200
-   :lineno-start: 190
+   :lines: 200-202
+   :lineno-start: 200
    :language: python
 
-The lines 203-262 produce the following pictures:
+The lines 205-2624produce the following pictures:
 
 
 .. image:: ../examples/MnSi/reference/TotalFields.png
@@ -160,10 +162,36 @@ The lines 203-262 produce the following pictures:
 
 Interestingly, left-handed and right-handed orders produce different local field. This is due to the lack of inversion symmetry.
 
-The phase
----------
+The phase between Mn atoms in the unit cell
+-------------------------------------------
+
+Dalmas De and Reotier colleagues pointed out that the complete description
+of the spiral phase requires three additional degrees of freedom: the two angles 
+describing the misalignmente between the rotation plane of the Mn moments
+and the plane perpendicular to the propagation vector :math:`k` and the
+phase between the helices in the Mn sites belonging to the two crystallographic 
+orbits (see [DalmasdeReotier2016]_).
+This last parameter strongly influence the results.
+
+Here we define a phase shift :math:`\phi=2` degrees as reported in the article 
+mentioned above.
+
+
+.. literalinclude:: ../examples/MnSi/run_example.py
+   :lines: 293-231
+   :lineno-start: 293
+   :emphasize-lines: 1,11
+   :language: python
+
+Repeating the same procedure discussed above leads to the following 
+comparison between the magnetic structures with :math:`\phi=0` (blue) and 
+and :math:`\phi=2` deg (orange).
+
+.. image:: ../examples/MnSi/reference/Histogram-zoom.png
+   :height: 150
+   :width: 300
+   :alt: Histogram of local fields
 
-TODO
 
 Bibliography
 ------------

examples/MnSi/run_example.py

@@ -18,9 +18,10 @@
 
 The magnetic structure of MnSi is characterized by spins forming a 
 left-handed incommensurate helix with a propagation vector k≃0.036 Å^−1
-in the [111] direction [5–7]. The static Mn moments ( ∼0.4μB for T→0 K)
+in the [111] direction [5–7]. The static Mn moments ( ∼0.385μB for T→0 K)
 point in a plane perpendicular to the propagation vector. 
 
+All data is taken from PRB 93 144419 (2016)
 """
 
 print("Create sample...", end='')
@@ -123,7 +124,7 @@
 a = 4.558 # Ang
 a_star = 2*np.pi/a #Ang^-1
 
-norm_k = 0.036 # Å −1
+norm_k = 0.035 # Å −1
 
 k_astar = k_bstar = k_cstar = (1/np.sqrt(3))*norm_k
 
@@ -145,10 +146,11 @@ def uv(vec):
                     # There are two choices here: left handed or right 
                     # handed spiral. We will do both.
 
-
+# Fianally, the experimental contact coupling term is:
+ContatExp = -0.066679616
 
 # assuming that, at x=0, a Mn moment is // a, the spiral can be obtained as
-RH_FCs = 0.4 * np.array([(a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)),
+RH_FCs = 0.385 * np.array([(a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)),
                          (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn2)),
                          (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn3)),
                          (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn4)),
@@ -158,7 +160,7 @@ def uv(vec):
                          [0,0,0]
                         ])
 # Notice the minus sign.
-LH_FCs = 0.4 * np.array([(a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)),
+LH_FCs = 0.385 * np.array([(a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)),
                          (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn2)),
                          (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn3)),
                          (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn4)),
@@ -196,8 +198,8 @@ def uv(vec):
 
 
 for i in range(4):
-    r_RH[i].ACont = -0.066
-    r_LH[i].ACont = -0.066
+    r_RH[i].ACont = ContatExp
+    r_LH[i].ACont = ContatExp
 
 
 # Two subplots, unpack the axes array immediately
@@ -227,8 +229,8 @@ def uv(vec):
 
 
 for i in range(4):
-    r_RH[i].ACont = -0.066
-    r_LH[i].ACont = -0.066
+    r_RH[i].ACont = ContatExp
+    r_LH[i].ACont = ContatExp
 
 # Two subplots, unpack the axes array immediately
 f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10,5))
@@ -247,6 +249,7 @@ def uv(vec):
     hist, bin_range = np.histogram(np.linalg.norm(r_LH[i].T, axis=1), bins=N_BINS, range=(0.08,0.24))
     LH_Hist += hist
 
+# just for plotting, gnerate intermediate positions for bins of the histogram
 mid_of_bin = bin_range[0:-1]+0.5*np.diff(bin_range)
 
 ax1.plot(mid_of_bin, RH_Hist)
@@ -260,3 +263,89 @@ def uv(vec):
 ax2.set_xlabel('B[T]')
 
 plt.savefig("Histogram.png")
+
+
+
+# Now add a phase between the two orbital types, see PRB 93 144419 (2016)
+
+# Two subplots, unpack the axes array immediately
+f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10,5))
+
+ax1.set_title('right-handed')
+ax2.set_title('left-handed')
+
+# replot old results fro perfect field to compare later
+ax1.plot(mid_of_bin, RH_Hist, label = 'original')
+ax2.plot(mid_of_bin, LH_Hist, label = 'original')
+
+ax1.set_ylim([0,1600])
+ax2.set_ylim([0,1600])
+ax1.set_xlim([0.08,0.1])
+ax2.set_xlim([0.08,0.1])
+
+ax1.set_ylabel('P(B)')
+ax1.set_xlabel('B[T]')
+ax2.set_xlabel('B[T]')
+
+
+
+# assuming that, at x=0, a Mn moment is // a, the new spiral can be obtained as
+RH_FCs = 0.385 * np.array([(a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)-np.pi*2./180.),
+                         (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn2)),
+                         (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn3)),
+                         (a+1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn4)),
+                         [0,0,0],
+                         [0,0,0],
+                         [0,0,0],
+                         [0,0,0]
+                        ])
+# Notice the minus sign.
+LH_FCs = 0.385 * np.array([(a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn1)-np.pi*2./180.),
+                         (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn2)),
+                         (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn3)),
+                         (a-1j*b)*np.exp(-2*np.pi*1j*np.dot(k_rlu,pos_Mn4)),
+                         [0,0,0],
+                         [0,0,0],
+                         [0,0,0],
+                         [0,0,0]
+                        ])
+
+s.new_mm()
+s.mm.desc = "Right handed spiral with 2 deg. pahse"
+s.mm.k = k_rlu
+s.mm.fc = RH_FCs
+
+s.new_mm()
+s.mm.desc = "Left handed spiral with 2 deg. pahse"
+s.mm.k = k_rlu
+s.mm.fc = LH_FCs
+
+
+s.current_mm_idx = 3;      # N.B.: indexes start from 0 but idx=0 is the transverse field!
+                           #       and 1 and 2 are the perfect spiral (no phase shift)
+r_RH = locfield(s, 'i',[50,50,50],100,nnn=3,nangles=36000)
+
+s.current_mm_idx = 4;
+r_LH = locfield(s, 'i',[50,50,50],100,nnn=3,nangles=36000)
+
+
+for i in range(4):
+    r_RH[i].ACont = ContatExp
+    r_LH[i].ACont = ContatExp
+
+
+LH_Hist=np.zeros(N_BINS)
+RH_Hist=np.zeros(N_BINS)
+for i in range(4):
+    hist, bin_range = np.histogram(np.linalg.norm(r_RH[i].T, axis=1), bins=N_BINS, range=(0.08,0.24))
+    RH_Hist += hist
+    hist, bin_range = np.histogram(np.linalg.norm(r_LH[i].T, axis=1), bins=N_BINS, range=(0.08,0.24))
+    LH_Hist += hist
+
+
+ax1.plot(mid_of_bin, RH_Hist, label='2 deg. shift')
+ax2.plot(mid_of_bin, LH_Hist, label='2 deg. shift')
+
+plt.legend()
+
+plt.savefig("Histogram-zoom.png")