Skip to content

Commit 5fea09a

Browse files
fangohrfangohr
committed
Merge branch 'master' of github.com:computationalmodelling/fidimag
2 parents dddc72d + 0bd9fc7 commit 5fea09a

File tree

9 files changed

+560
-103
lines changed

9 files changed

+560
-103
lines changed

‎examples/micromagnetic/dw_stt_cpp/main.py‎

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
import matplotlib as mpl
2+
mpl.use("Agg")
3+
import matplotlib.pyplot as plt
4+
import matplotlib.gridspec as gridspec
5+
6+
import numpy as np
7+
from fidimag.micro import Sim, UniformExchange, UniaxialAnisotropy
8+
from fidimag.common import CuboidMesh
9+
import fidimag.common.constant as const
10+
11+
from scipy.optimize import curve_fit
12+
13+
14+
mu0 = 4 * np.pi * 1e-7
15+
global_Ms = 8.6e5
16+
global_A = 1.3e-11
17+
global_K = 0.25*mu0*global_Ms**2
18+
19+
global_P = 0.5
20+
global_d = 2e-9 #2nm
21+
global_const = const.h_bar*const.gamma*global_P/(2*const.c_e*global_d*global_Ms)
22+
23+
def init_m(pos):
24+
25+
x = pos[0]
26+
27+
if x < 300:
28+
return (1, 0, 0)
29+
elif 300 <= x < 400:
30+
return (0, 1, 1)
31+
else:
32+
return (-1, 0, 0)
33+
34+
35+
def relax_system(mesh):
36+
37+
# Only relaxation
38+
sim = Sim(mesh, name='relax')
39+
40+
# Simulation parameters
41+
sim.set_tols(rtol=1e-8, atol=1e-10)
42+
sim.alpha = 0.5
43+
sim.gamma = 2.211e5
44+
sim.Ms = 8.6e5
45+
sim.do_precession = False
46+
47+
# The initial state passed as a function
48+
sim.set_m(init_m)
49+
# sim.set_m(np.load('m0.npy'))
50+
51+
# Energies
52+
A = 1.3e-11
53+
exch = UniformExchange(A=A)
54+
sim.add(exch)
55+
56+
anis = UniaxialAnisotropy(5e4)
57+
sim.add(anis)
58+
59+
# Start relaxation and save the state in m0.npy
60+
sim.relax(dt=1e-14, stopping_dmdt=0.00001, max_steps=5000,
61+
save_m_steps=None, save_vtk_steps=None)
62+
63+
np.save('m0.npy', sim.spin)
64+
65+
66+
67+
def model(xs, x0):
68+
delta = 16.1245154966
69+
return -np.tanh((xs-x0)/delta)
70+
71+
def extract_dw(m):
72+
73+
m.shape=(-1,3)
74+
mx = m[:,0]
75+
my = m[:,1]
76+
mz = m[:,2]
77+
78+
xs = np.linspace(0.5, 999.5, 1000)
79+
80+
id = np.argmin(abs(mx))
81+
82+
popt, pcov = curve_fit(model, xs, mx, p0=[id*1.0])
83+
print(popt[0], id)
84+
85+
theta = np.arctan2(mz[id], my[id])
86+
return popt[0], theta
87+
88+
def sinc_fun(t):
89+
return np.sinc(1e9*t) # sinc(x) = sin(\pi x)/(\pi x)
90+
91+
# This function excites the system with a
92+
# current in the x-direction using Spin Transfer Torque
93+
# formalism
94+
def excite_system(mesh, beta=0.0):
95+
96+
# Specify the stt dynamics in the simulation
97+
sim = Sim(mesh, name='dyn_%g'%beta, driver='llg_stt_cpp')
98+
99+
sim.set_tols(rtol=1e-12, atol=1e-12)
100+
sim.alpha = 0.1
101+
sim.gamma = 2.211e5
102+
sim.Ms = 8.6e5
103+
104+
# sim.set_m(init_m)
105+
sim.set_m(np.load('m0.npy'))
106+
107+
# Energies
108+
A = 1.3e-11
109+
exch = UniformExchange(A=A)
110+
sim.add(exch)
111+
112+
anis = UniaxialAnisotropy(5e4)
113+
sim.add(anis)
114+
115+
# beta is the parameter in the STT torque
116+
sim.a_J = global_const*1e11
117+
sim.p = (1,0,0)
118+
sim.beta = beta
119+
120+
# The simulation will run for 5 ns and save
121+
# 500 snapshots of the system in the process
122+
ts = np.linspace(0, 0.5e-9, 21)
123+
124+
xs=[]
125+
thetas=[]
126+
127+
for t in ts:
128+
print('time', t)
129+
sim.run_until(t)
130+
spin = sim.spin.copy()
131+
x, theta = extract_dw(spin)
132+
xs.append(x)
133+
thetas.append(theta)
134+
sim.save_vtk()
135+
136+
np.savetxt('dw_%g.txt'%beta,np.transpose(np.array([ts, xs,thetas])))
137+
138+
def analytical(ts):
139+
delta = 16.1245154966
140+
alpha = 0.1
141+
beta = 0.16
142+
aj = global_const*1e11
143+
theta = (beta-alpha)/(1+alpha**2)*ts*aj
144+
z = (1+alpha*beta)/(1+alpha**2)*delta*aj*ts
145+
return z, theta
146+
147+
148+
def plot():
149+
data=np.loadtxt('dw_0.16.txt')
150+
ts = data[:,0]
151+
xs = data[:,1]-data[0,1]
152+
theta = data[:,2]-data[0,2]
153+
154+
z_an, theta_an = analytical(ts)
155+
156+
fig = plt.figure(figsize=(5, 2.6))
157+
158+
gs = gridspec.GridSpec(1, 2, width_ratios=[4, 4])
159+
ax0 = plt.subplot(gs[0])
160+
ax0.plot(ts*1e9,xs,'.')
161+
ax0.plot(ts*1e9, z_an, '-')
162+
plt.ylabel(r'DW shift')
163+
plt.xlabel(r'Time (ns)')
164+
165+
166+
ax1 = plt.subplot(gs[1])
167+
ax1.plot(ts*1e9,theta,'.')
168+
ax1.plot(ts*1e9, theta_an, '-')
169+
plt.ylabel(r'$\phi$')
170+
plt.xlabel(r'Time (ns)')
171+
172+
plt.tight_layout()
173+
analytical(ts)
174+
175+
plt.savefig('xs_theta.pdf')
176+
177+
178+
if __name__ == '__main__':
179+
180+
# We will crate a mesh with 1000 elements of elements
181+
# in the x direction, and 1 along y and z
182+
# (so we have a 1D system)
183+
mesh = CuboidMesh(nx=1000, ny=1, nz=1,
184+
dx=1.0, dy=2, dz=2.0,
185+
unit_length=1e-9)
186+
187+
# Relax the initial state
188+
relax_system(mesh)
189+
190+
excite_system(mesh, beta=0.16)
191+
plot()
192+

‎fidimag/atomistic/anisotropy.py‎

Lines changed: 36 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,18 +7,49 @@
77
class Anisotropy(Energy):
88

99
"""
10-
compute the anisotropy field with the energy density E = - K (m.u)^2
10+
11+
This class provides an anisotropy term to the system energy, which
12+
is defined as
13+
14+
__ -> ^ 2
15+
E = - \ K_i ( S_i * u_i )
16+
/__
17+
i
18+
19+
with K_i the anisotropy constant at the i-th site and u_i the unitary
20+
direction of the anisotropy vector at the i-th site, thus the magnitude and
21+
axis can be space dependent.
22+
23+
OPTIONAL ARGUMENTS --------------------------------------------------------
24+
25+
Ku :: The anisotropy constant. Can be a constant or a space
26+
dependent scalar field, given as a function or array.
27+
28+
axis :: The unitary axis vector. It can be a 3 tuple, list
29+
or array (uniaxial anisotropy), or a space dependent
30+
vector field, gicen as a function or array.
31+
32+
name :: Interaction name
33+
34+
USAGE ---------------------------------------------------------------------
35+
36+
Considering a simulation object *Sim*, an uniaxial anisotropy along
37+
the z direction can be defined as
38+
39+
K = 1 * meV
40+
Sim.add(Anisotropy(K, axis=(0, 0, 1)))
41+
42+
For a space dependent anisotropy along the x direction, we can define a
43+
scalar field for the magnitude and a vector field for the anisotropy axes.
44+
1145
"""
1246

13-
def __init__(self, Ku, axis=(1, 0, 0), name='anis', direction=None):
47+
def __init__(self, Ku, axis=(1, 0, 0), name='anis'):
1448
self.Ku = Ku
1549
self.name = name
1650
self.axis = axis
1751
self.jac = True
1852

19-
if direction is not None:
20-
self.axis = direction
21-
2253
def setup(self, mesh, spin, mu_s):
2354
super(Anisotropy, self).setup(mesh, spin, mu_s)
2455

‎fidimag/atomistic/dmi.py‎

Lines changed: 72 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -6,11 +6,74 @@
66
class DMI(Energy):
77

88
"""
9-
Hamiltonian = D*[S_i x S_j]
10-
==> H = D x S_j
9+
10+
This class provides the Dzyaloshinskii-Moriya Interaction (DMI) energy
11+
term, defined as
12+
13+
__ -> -> ->
14+
E = \ D_ij * ( S_i X S_j )
15+
/__
16+
<i, j>
17+
i != j
18+
19+
where D_ij is the Dzyaloshinskii vector, S_i and S_j are the total spin
20+
vectors at the i-th and j-th lattice sites, and <i, j> means counting the
21+
interaction between neighbouring spins only once. The Dzyaloshinskii vector
22+
magnitude can vary with space.
23+
24+
Currently, there are implemented two kinds of DMI, that depend on the
25+
Dzyaloshinskii vector definition. Calling D_i the Dzyaloshinskii vector
26+
magnitude at the i-th site:
27+
28+
BULK:
29+
-> ^
30+
D_ij = D_i r_ij with r_ij as the unit vector pointing from the i-th
31+
towards the j-th site
32+
(only working for square lattices)
33+
34+
INTERFACIAL:
35+
36+
-> ^ ^
37+
D_ij = D_i ( r_ij X z ) with r_ij as the unit vector pointing from the
38+
i-th towards the j-th site and z the unitary
39+
vector perpendicular to the magnetic material
40+
plane which, in theory, is in contact with a
41+
material with larger SOC (e.g. a metal).
42+
Thus, this DMI is defined in 2D, for
43+
interfaces or thin films.
44+
45+
For further details about the DMI calculations, take a look
46+
to the C library documentation (dmi.c)
47+
48+
49+
OPTIONAL ARGUMENTS: -------------------------------------------------------
50+
51+
dmi_type :: 'bulk' or 'interfacial'
52+
name :: Interaction name
53+
54+
USAGE: --------------------------------------------------------------------
55+
56+
If the DMI is homogeneous, it can be specified in a simulation object
57+
*Sim* as
58+
59+
Sim.add(DMI(D, dmi_type='bulk'))
60+
61+
where D is a float.
62+
63+
Otherwise, it can be specified as any Fidimag scalar field, passing a
64+
function or an array. For example, a space dependent DMI that changes
65+
linearly in the x-direction can be defined as:
66+
67+
def my_DMI(pos):
68+
D = 0.01 * meV
69+
return D * pos[0]
70+
71+
# Add DMI to Simulation object
72+
Sim.add(DMI(my_DMI))
73+
1174
"""
1275

13-
def __init__(self, D, name='dmi', dmi_type='bulk'):
76+
def __init__(self, D, name='DMI', dmi_type='bulk'):
1477
self.D = D
1578

1679
self.name = name
@@ -32,20 +95,20 @@ def setup(self, mesh, spin, mu_s):
3295
# We will generate the Dzyaloshinskii vectors according
3396
# to the lattice, for the Interfacial DMI
3497
self.DMI_vector = self.compute_DMI_vectors(self.nneighbours)
35-
98+
3699
if self.dmi_type == 'bulk':
37100
self._D = np.zeros(self.neighbours.shape)
38101
if isinstance(self.D, (int, float)):
39-
self._D[:,:] = self.D
102+
self._D[:, :] = self.D
40103
elif hasattr(self.D, '__call__'):
41-
n = self.mesh.n
104+
n = self.mesh.n
42105
for i in range(n):
43106
value = self.D(self.coordinates[i])
44107
if len(value) == 6:
45-
self._D[i,:] = value[:]
108+
self._D[i, :] = value[:]
46109
else:
47110
raise Exception('The given spatial function for D is not acceptable!')
48-
111+
49112

50113
def compute_field(self, t=0, spin=None):
51114

@@ -55,7 +118,7 @@ def compute_field(self, t=0, spin=None):
55118
m = self.spin
56119

57120
if self.dmi_type == 'bulk':
58-
clib.compute_dmi_field(m,
121+
clib.compute_dmi_field(m,
59122
self.field,
60123
self.energy,
61124
self._D,

0 commit comments

Comments
 (0)