@@ -87,32 +87,51 @@ void dmi_field_bulk(double *m, double *field, double *energy, double *Ms_inv,
8787 * where j is the index of a NN neighbour and D_{ij} = D * r_{ij}, r_{ij}
8888 * being the position vector connectin the site i with the site j
8989 *
90- * In the continuum, the field can be written as : - nabla X M , which can
91- * be discretised with an expression similar to the atomic model one when
92- * doing the finite difference approximation of the derivative, thus we
93- * only need the DMI vector as: -D * r_{ij}
94- *
95- * So, in the loop through neighbours we compute:
96- *
97- * neighbour field sum this gives
98- * -x: (D / dx) * (+x X M) --> Components in y, z
99- * +x: (D / dx) * (-x X M) --> %
100- * -y: (D / dy) * (+y X M) --> Components in x, z
101- * +y: (D / dy) * (-y X M) --> %
102- * -z: (D / dz) * (+z X M) --> Components in x, y
103- * +z: (D / dz) * (-z X M) --> %
104- *
105- * which gives:
106- * field_x = (D / dx) * (mz[-x] - mz[+x])
107- * + (D / dz) * (mx[-z] - mx[+x])
90+ * In the continuum, the field can be written as (M = Ms m):
91+ * ->
92+ * - (2 D / mu0 Ms) nabla X m ,
93+ *
94+ * which can be discretised with an expression similar to the atomic model
95+ * one when doing the finite difference approximation of the derivative,
96+ * thus we only need the DMI vector as: -D * r_{ij}. Specifically, we have
97+ * that:
98+ *
99+ * ----> -> ^ ^
100+ * field = - 2 D nabla X m = -D [ (d m_z - d m_y) x + (...) y + ]
101+ * ---- ----- -----
102+ * mu0 Ms d y d z
103+ *
104+ * Discretising the derivatives, the field components are:
105+ *
106+ * field_x = - (2 D / mu0 Ms) [ (1 / 2 dy) * (mz[+y] - mz[-y]) -
107+ * (1 / 2 dz) * (my[+z] - my[-z])
108+ * ]
109+ *
110+ * field_y = - (2 D / mu0 Ms) [ (1 / 2 dz) * (mx[+z] - mx[-z]) -
111+ * (1 / 2 dx) * (mz[+x] - mz[-x])
112+ * ]
108113 * ...
109114 *
110- * which are the first order derivatives.
115+ * where +x, -x, etc are the neighbours at position +x, -x, etc
116+ * We now can collect the terms involving the neighbours at
117+ * +x, -x, +y, -y, +z and -z.
118+ * For example the field H where the +x neighbour is involved is:
119+ *
120+ * -> ^ ^
121+ * H(+x) = (-2 D / mu0 Ms) * ( my(+x) z - mz(+x) y ) * (1 / 2 dx)
122+ *
123+ * ^ ->
124+ * = (- D / mu0 Ms dx) * ( x X m )
125+ *
126+ * In general, the field for every component will be the cross product of
127+ * r_ij with m, where r_ij is the vector towards the j-th neighbour, as in
128+ * the discrete case. This term is divided by the mesh discretisation
129+ * (dx, dy or dz) in the corresponding direction
111130 *
112131 * The DMI vector norms are given by the *D array. If our simulation has
113132 * n mesh nodes, then the D array is (6 * n) long, i.e. the DMI vector norm
114133 * per every neighbour per every mesh node. The order is the same than the NNs
115- * array, i.e.
134+ * array, i.e.
116135 *
117136 * D = [D(x)_0, D(-x)_0, D(y)_0, D(-y)_0, D(z)_0, D(-z)_0, D(x)_1, ...]
118137 *
@@ -130,13 +149,14 @@ void dmi_field_bulk(double *m, double *field, double *energy, double *Ms_inv,
130149 /* The DMI vector directions are the same according to the neighbours
131150 * positions. Thus we set them here to avoid compute them
132151 * every time in the loop . So, if we have the j-th NN,
133- * the DMI vector will be dmivector[3 * j] */
134- double dmivector [18 ] = {-1 , 0 , 0 ,
135- 1 , 0 , 0 ,
136- 0 , -1 , 0 ,
137- 0 , 1 , 0 ,
138- 0 , 0 , -1 ,
139- 0 , 0 , 1
152+ * the DMI vector will be dmivector[3 * j]
153+ * */
154+ double dmivector [18 ] = {-1 , 0 , 0 , // -x
155+ 1 , 0 , 0 , // +x
156+ 0 , -1 , 0 , // -y
157+ 0 , 1 , 0 , // +y
158+ 0 , 0 , -1 , // -z
159+ 0 , 0 , 1 // +z
140160 };
141161
142162 /* These are for the DMI prefactor or coefficient */
@@ -175,7 +195,7 @@ void dmi_field_bulk(double *m, double *field, double *energy, double *Ms_inv,
175195 * is larger than zero */
176196 if (Ms_inv [ngbs [idn + j ]] > 0 ){
177197
178- /* We do here: (D / dx_i) * ( r_{ij} X M_{j} )
198+ /* We do here: - (D / dx_i) * ( r_{ij} X M_{j} )
179199 * The cross_i function gives the i component of the
180200 * cross product. The coefficient is computed according
181201 * to the DMI strength of the current lattice site.
@@ -232,44 +252,64 @@ void dmi_field_interfacial(double *m, double *field, double *energy, double *Ms_
232252 /* In the atomic model, the effective field in the i-th spin site has the
233253 * summation: D_{ij} x S_{ij}
234254 *
235- * For the Interfacial DMI, D_{ij} has the structure: D_{ij} = r_{ij} X z
236- * See Rohart et al. Phys. Rev. B 88, 184422)
237- * We will use R&T's convention for the DMI sign.
255+ * For the Interfacial DMI, we use the following energy density structure
256+ * with Lifshitz invariants (see Rohart et al. Phys. Rev. B 88, 184422))
257+ *
258+ * w_DM = D ( L_{xz}^{(x)} + L_{yz}^{(y)} )
259+ *
260+ * Hence, using a variational derivative, the field can be computed as
261+ *
262+ * / -> -> \
263+ * field = - 2 D | ^ dm ^ dm |
264+ * ---- | y X -- - x X -- |
265+ * mu0 Ms \ dx dy /
266+ *
267+ * Using finite differences for the derivatives, we can get the components
268+ * contributed to the field by the neighbours in the +x, -x, +y and -y
269+ * directions. For instance,
238270 *
239- * Notice that in [Yang et al. Phys. Rev. Lett. 115, 267210] they use the
240- * opposite sign.
271+ * / \
272+ * field(+x) = -2 D | mz(+x) ^ mx(+x) ^ |
273+ * ---- | ----- x - ----- z |
274+ * mu0 Ms \ 2 dx 2 dx /
275+ *
276+ * - D 1 ^ ->
277+ * = ---- -- ( y X M(+x) )
278+ * mu0 Ms dx
241279 *
242- * The Dzyaloshinskii vectors are IN plane. This function only works
243- * along the XY plane since we assume there is a non magnetic material below
244- * with a different SOC which gives the DMI for neighbouring in plane spins
245280 *
246- * The computation is similar than before, only that we no longer have
247- * a z component. In the continuum the H_dmi field is:
281+ * In general, the derivative will be given as in the discrete spin model
282+ * by using the Dzyaloshinskii vector as
248283 *
249- * H_dmi = (2 * D / (Ms a a_z)) * ( x X dM/dy - y X dM/dx )
284+ * ^ ->
285+ * D_ij = z X r_ij
250286 *
251- * where "a" is the lattice spacing in the plane and "a_z" the system
252- * thickness along z
287+ * and dividing by the mesh discretisation in the direction of the
288+ * neighbour.
289+ * This DMI vector convention is given by [Yang et al. PRL 115, 267210]
253290 *
254- * This derivative can be obtained doing the cross product
255- * as in the atomic model, hence when going through the neighbours
291+ * The Dzyaloshinskii vectors are IN plane. This function only works along
292+ * the XY plane since we assume there is a non magnetic material below with
293+ * a different SOC which gives the DMI for neighbouring in plane spins
294+ *
295+ * As in the atomic model, when going through the neighbours
256296 * loop, we can obtain the field doing the following cross products:
257297 *
258298 * neighbour field sum
259- * -x: (D / dx) * (+ y X M)
260- * +x: (D / dx) * (- y X M)
261- * -y: (D / dy) * (- x X M)
262- * +y (D / dy) * (+ x X M)
299+ * -x: (D / dx) * (- y X M)
300+ * +x: (D / dx) * (+ y X M)
301+ * -y: (D / dy) * (+ x X M)
302+ * +y (D / dy) * (- x X M)
263303 *
264304 * So, our Dzyaloshinskii vectors can be depicted in the cuboid mesh as
265305 *
266306 * o +y
267307 * |
268- * --> D
269- * ^ |
308+ * <-- D
309+ * | ^
270310 * -x o __ | __ o __ | _ o +x
271- * | v
272- * <--
311+ * v |
312+ * -->
273313 * |
274314 * o -y
275315 *
@@ -279,7 +319,7 @@ void dmi_field_interfacial(double *m, double *field, double *energy, double *Ms_
279319 * The DMI vector norms are given by the *D array. If our simulation has
280320 * n mesh nodes, then the D array is (4 * n) long, i.e. the DMI vector norm
281321 * per every neighbour per every mesh node. The order is the same than the NNs
282- * array, i.e.
322+ * array, i.e.
283323 *
284324 * D = [D(x)_0, D(-x)_0, D(y)_0, D(-y)_0, D(x)_1, ...]
285325 *
@@ -294,10 +334,10 @@ void dmi_field_interfacial(double *m, double *field, double *energy, double *Ms_
294334 * (NNs are in the order: -x, +x, -y, +y)
295335 * For interfacial DMI we only compute the DMI in 2D
296336 * */
297- double dmivector [12 ] = { 0 , 1 , 0 ,
298- 0 , -1 , 0 ,
299- -1 , 0 , 0 ,
337+ double dmivector [12 ] = { 0 , -1 , 0 ,
338+ 0 , 1 , 0 ,
300339 1 , 0 , 0 ,
340+ -1 , 0 , 0 ,
301341 };
302342 double dxs [4 ] = {dx , dx , dy , dy };
303343
@@ -327,7 +367,7 @@ void dmi_field_interfacial(double *m, double *field, double *energy, double *Ms_
327367 if (ngbs [idn + j ] >= 0 ) {
328368
329369 /* DMI coefficient according to the neighbour position */
330- DMIc = D [idn_DMI + j ] / dxs [j ];
370+ DMIc = - D [idn_DMI + j ] / dxs [j ];
331371
332372 /* Magnetisation array index of the neighbouring spin
333373 * since ngbs gives the neighbour's index */
0 commit comments