diff --git a/cmtj/models/sb.py b/cmtj/models/sb.py index 4dadbbf..f48c5c1 100644 --- a/cmtj/models/sb.py +++ b/cmtj/models/sb.py @@ -174,6 +174,49 @@ def __any_grad_iec_interaction(self, J, layer: "LayerSB"): return [dEdtheta, dEdphi, d2Edtheta2, d2Edphi2, d2Edphidtheta] + def dmi_interaction(self, D, layer: "LayerSB"): + """ + DMi energy is Edmi = D z(m1 x m2) + """ + if layer is None: + return 0 + [m1x, m1y, _] = self.m.get_cartesian() + [m2x, m2y, _] = layer.m.get_cartesian() + return D * (m1x * m2y - m1y * m2x) + + def grad_dmi_interaction(self, D_top, D_bottom, top_layer: "LayerSB", + bottom_layer: "LayerSB"): + """ + DMI energy and its gradient are both not symmetric and is computed per layer + Compute joint interaction on the layers from the top and bottom + """ + grad1 = self.__any_grad_dmi_interaction(D_top, top_layer) + grad2 = self.__any_grad_dmi_interaction(D_bottom, bottom_layer) + return [ + grad1[0] + grad2[0], grad1[1] + grad2[1], grad1[2] + grad2[2], + grad1[3] + grad2[3], grad1[4] + grad2[4] + ] + + def __any_grad_dmi_interaction(self, D, layer: "LayerSB"): + if (layer is None) or D == 0: + return [0, 0, 0, 0, 0] + constD = D + + dEdtheta = -constD * math.sin( + layer.m.theta) * math.sin(self.m.phi - layer.m.phi) * self.ctheta + + dEdphi = -constD * self.stheta * math.sin( + layer.m.theta) * math.cos(self.m.phi - layer.m.phi) + + d2Edtheta2 = constD * self.stheta * math.sin(self.m.phi - layer.m.phi) + + d2Edphi2 = d2Edtheta2 + + d2Edphidtheta = -D * math.sin( + layer.m.phi) * self.ctheta * math.cos(self.m.phi - layer.m.phi) + + return [dEdtheta, dEdphi, d2Edtheta2, d2Edphi2, d2Edphidtheta] + def surface_anisotropy(self): return (-self.Ks + 0.5 * self.Ms * mu0) * self.ctheta**2 @@ -213,22 +256,26 @@ def grad_volume_anisotropy(self): return [dEdtheta, dEdphi, d2Edtheta2, d2Edphi2, d2Edphidtheta] def compute_grad_energy(self, Hinplane: VectorObj, Jtop: float, - Jbottom: float, top_layer: "LayerSB", - bottom_layer: "LayerSB"): + Jbottom: float, Dtop: float, Dbottom: float, + top_layer: "LayerSB", bottom_layer: "LayerSB"): g1 = self.grad_ext_field(Hinplane) g2 = self.grad_surface_anisotropy() g3 = self.grad_volume_anisotropy() g4 = self.grad_iec_interaction(Jtop, Jbottom, top_layer, bottom_layer) - return [g1[i] + g2[i] + g3[i] + g4[i] for i in range(5)] + g5 = self.grad_dmi_interaction(Dtop, Dbottom, top_layer, bottom_layer) + return [g1[i] + g2[i] + g3[i] + g4[i] + g5[i] for i in range(5)] def compute_energy(self, Hinplane: VectorObj, Jtop: float, Jbottom: float, - top_layer: "LayerSB", bottom_layer: "LayerSB"): + Dtop: float, Dbottom: float, top_layer: "LayerSB", + bottom_layer: "LayerSB"): e1 = self.ext_field(Hinplane) e2 = self.surface_anisotropy() e3 = self.volume_anisotropy() e4 = self.iec_interaction( Jbottom, bottom_layer) + self.iec_interaction(Jtop, top_layer) - evec = [e1, e2, e3, e4] + e5 = self.dmi_interaction(Dtop, top_layer) + self.dmi_interaction( + Dbottom, bottom_layer) + evec = [e1, e2, e3, e4, e5] return evec def get_current_position(self): @@ -239,7 +286,8 @@ def set_current_position(self, pos): self.m.phi = pos[1] def compute_frequency_at_equilibrum(self, Hinplane: VectorObj, Jtop: float, - Jbottom: float, top_layer: "LayerSB", + Jbottom: float, Dtop: float, + Dbottom: float, top_layer: "LayerSB", bottom_layer: "LayerSB"): """Computes the resonance frequency (FMR) of the layers :param Hinplance: vector that describes the applied H. @@ -248,8 +296,9 @@ def compute_frequency_at_equilibrum(self, Hinplane: VectorObj, Jtop: float, :param top_layer: LayerSB definition of the layer above the current one. :param bottom layer: LayerSB definition of the layer below the current one.""" (_, _, d2Edtheta2, d2Edphi2, - d2Edphidtheta) = self.compute_grad_energy(Hinplane, Jbottom, Jtop, - top_layer, bottom_layer) + d2Edphidtheta) = self.compute_grad_energy(Hinplane, Jtop, Jbottom, + Dtop, Dbottom, top_layer, + bottom_layer) if self.stheta != 0.: fmr = (d2Edtheta2 * d2Edphi2 - d2Edphidtheta**2) / math.pow( @@ -271,16 +320,19 @@ def __init__(self, layers: List[LayerSB], Hext: VectorObj, J: List[float] = [], + D: List[float] = [], silent: bool = False) -> None: """ :param layers: list of LayerSB, layer definitions. :param Hext: applied external magnetic field vector. Defined for all the layers. - :param J: the list of IEC constants. 0 element defines coupling between 0 and 1 layer, etc. + :param J: list of IEC constants. 0 element defines coupling between 0 and 1 layer, etc. + :param D: list of DMI constants. :param silent: debug mode? defaults to False. """ if len(layers) != (len(J) + 1): raise ValueError("Number of layers must be equal to len(J) + 1") self.J = J + self.D = D self.layers: List[LayerSB] = layers self.energy = np.zeros((len(self.layers), ), dtype=np.float32) self.current_position = np.zeros((len(self.layers), 2), @@ -291,30 +343,41 @@ def __init__(self, self.ener_labels = ['external', 'surface', 'volume', 'iec'] self.history = defaultdict(list) + def __get_interaction_constant(self, layers: List[LayerSB], + layer_indx: int, + interaction_constant_list: List[float]): + """Simple function to figure out top and bottom handles and interaction constants""" + if layer_indx > 0: + bottom_handle = layers[layer_indx - layer_indx] + interaction_const_bottom = interaction_constant_list[layer_indx - + 1] + else: + bottom_handle = None + interaction_const_bottom = 0 + if layer_indx < len(layers) - 1: + top_handle = layers[layer_indx + 1] + interaction_const_top = interaction_constant_list[layer_indx] + else: + top_handle = None + interaction_const_top = 0 + return interaction_const_top, interaction_const_bottom, top_handle, bottom_handle + def compute_energy_step(self): """ Note that symmetric coupling is only computed once and not per layer """ for i, layer in enumerate(self.layers): - if i > 0: - bottom_layer_handle = self.layers[i - 1] - Jbottom = self.J[i - 1] - else: - bottom_layer_handle = None - Jbottom = 0 - if i < len(self.layers) - 1: - top_layer_handle = self.layers[i + 1] - Jtop = self.J[i] - else: - top_layer_handle = None - Jtop = 0 - - evec = layer.compute_energy(self.Hext, Jtop, Jbottom, - top_layer_handle, bottom_layer_handle) + Jtop, Jbottom, top_layer_handle, bottom_layer_handle = self.__get_interaction_constant( + self.layers, i, self.J) + Dtop, Dbottom, _, _ = self.__get_interaction_constant( + self.layers, i, self.D) + evec = layer.compute_energy(self.Hext, Jtop, Jbottom, Dtop, + Dbottom, top_layer_handle, + bottom_layer_handle) self.energy[i] = sum(evec) self.current_position[i, :] = layer.get_current_position() self.grad_energy[i, :] = layer.compute_grad_energy( - self.Hext, Jtop, Jbottom, top_layer_handle, + self.Hext, Jtop, Jbottom, Dtop, Dbottom, top_layer_handle, bottom_layer_handle) for ener_val, ener_label in zip(evec, self.ener_labels): self.history[f"energy_{ener_label}_{i}"].append(ener_val) @@ -353,20 +416,13 @@ def gradient_descent(self, def get_fmr(self, layer_index: int) -> float: """Computes FMR values for each of the layers in the model.""" - if layer_index > 0: - bottom_layer_handle = self.layers[layer_index - 1] - Jbottom = self.J[layer_index - 1] - else: - bottom_layer_handle = None - Jbottom = 0 - if layer_index < len(self.layers) - 1: - top_layer_handle = self.layers[layer_index + 1] - Jtop = self.J[layer_index] - else: - top_layer_handle = None - Jtop = 0 + Jtop, Jbottom, top_layer_handle, bottom_layer_handle = self.__get_interaction_constant( + self.layers, layer_index, self.J) + Dtop, Dbottom, _, _ = self.__get_interaction_constant( + self.layers, layer_index, self.D) return self.layers[layer_index].compute_frequency_at_equilibrum( - self.Hext, Jtop, Jbottom, top_layer_handle, bottom_layer_handle) + self.Hext, Jtop, Jbottom, Dtop, Dbottom, top_layer_handle, + bottom_layer_handle) def print_summary(self, steps): """Prints summary of the solution if silent param is true."""