|
| 1 | +from __future__ import print_function |
| 2 | +import pytest |
| 3 | + |
| 4 | +""" |
| 5 | +
|
| 6 | +Test of the NEBM on an elongated particle of cylindrical shape, based on [1], |
| 7 | +where a series of simple magnetic systems are tested using the NEBM in |
| 8 | +Spherical coordinates with an Euclidean distance. For the elongated particle |
| 9 | +test, the system is a 70 nm long and 12 nm wide (originally it is 13 nm) |
| 10 | +nanocylinder, with a uniaxial anisotropy in the direction of the long axis.The |
| 11 | +test goal is to find the minimum energy path (MEP) between two degenerate |
| 12 | +ferromagnetic states (saturated states along the anisotropy axis), which are |
| 13 | +the ground states of this system. The MEP is given by a transverse domain wall |
| 14 | +propagation that reverse the ferromagnetic state towards the opposite |
| 15 | +anisotropic direction. |
| 16 | +
|
| 17 | +In the test, we use the Geodesic and Cartesian NEBM codes (there are issues |
| 18 | +with the definition of the angles in the Spherical code). For the Geodesic code |
| 19 | +we also describe the magnetisation in Cartesian coordinates but we use a |
| 20 | +Geodesic distance as in [2] and a projection of the tangents. In the case of |
| 21 | +the Cartesian code we use an Euclidean distance and we do not use projections. |
| 22 | +The initial energy band we provide for the algorithm is similar to a coherent |
| 23 | +rotation of the spins, which is not the optimal path. |
| 24 | +
|
| 25 | +For this particular code, the cylinder long axis (and so the anisotropy axis) |
| 26 | +is defined along the z direction. In case we used spherical coordinates, |
| 27 | +the angles that define the spin directions are not completely defined when |
| 28 | +the spins point in the z direction (poles) thus the algorithm struggles |
| 29 | +to find the minimum energy path. It seems this is also the case for Cartesian |
| 30 | +coordinates since we haven't projected the tangents. |
| 31 | +
|
| 32 | +An interesting test would be to compare the convergences when using a more |
| 33 | +symmetrical initial band, like |
| 34 | + initial_images = [[0, 0, -1], [0, 1, 0], [0, 0, 1]] |
| 35 | +
|
| 36 | +According to our tests, the Geodesic code is the one that converges faster |
| 37 | +and not having much issues with the spin directions. |
| 38 | +
|
| 39 | +[1] Dittrich et al., JMMM 250 (2002) L12-L19 |
| 40 | +[2] Bessarab et al., Computer Physics Communications 196 (2015) 335-347 |
| 41 | +
|
| 42 | +""" |
| 43 | + |
| 44 | +# FIDIMAG: |
| 45 | +from fidimag.micro import Sim |
| 46 | +from fidimag.common import CuboidMesh |
| 47 | +from fidimag.micro import UniformExchange, UniaxialAnisotropy, Demag |
| 48 | +from fidimag.common.nebm_geodesic import NEBM_Geodesic |
| 49 | +from fidimag.common.nebm_cartesian import NEBM_Cartesian |
| 50 | +import numpy as np |
| 51 | + |
| 52 | +# Material Parameters --------------------------------------------------------- |
| 53 | + |
| 54 | +A = 10e-12 |
| 55 | +Kx = 3e5 |
| 56 | +Ms = 3.98e5 |
| 57 | + |
| 58 | +# ----------------------------------------------------------------------------- |
| 59 | + |
| 60 | + |
| 61 | +def relax_neb(sim, k, maxst, simname, initial_images, interpolations, |
| 62 | + save_every=10000, method='Geodesic', |
| 63 | + interpolation_method='linear' |
| 64 | + ): |
| 65 | + """ |
| 66 | + Relax a Fidimag simulation using the NEBM: |
| 67 | +
|
| 68 | + sim :: Simulation object with the system specifications |
| 69 | + k :: NEBM spring constant |
| 70 | + maxst :: NEBM max steps for the algorithm evolution |
| 71 | + simname :: NEBM simulation name |
| 72 | + initial_images :: NEBM images |
| 73 | + interpolations :: A list with the number of inteprolations between |
| 74 | + the NEBM images |
| 75 | + save_every :: Save VTK or NPY files every certain number of steps |
| 76 | + method :: NEBM coordinates: 'Cartesian', 'Geodesic' |
| 77 | + interpolation_method :: Method for interpolating the initial_images, |
| 78 | + 'linear' or 'rotation' (Rodrigues formulae) |
| 79 | + """ |
| 80 | + |
| 81 | + method_dict = {'Cartesian': NEBM_Cartesian, 'Geodesic': NEBM_Geodesic} |
| 82 | + |
| 83 | + neb = method_dict[method](sim, |
| 84 | + initial_images, |
| 85 | + interpolations=interpolations, |
| 86 | + spring_constant=k, |
| 87 | + name=simname, |
| 88 | + interpolation_method='rotation', |
| 89 | + openmp=True |
| 90 | + ) |
| 91 | + |
| 92 | + neb.relax(max_iterations=maxst, |
| 93 | + save_vtks_every=save_every, |
| 94 | + save_npys_every=save_every, |
| 95 | + stopping_dYdt=1) |
| 96 | + |
| 97 | + # x, E_interp = neb.compute_polynomial_approximation(200) |
| 98 | + # np.savetxt('interpolated_band_GEODESIC.dat', |
| 99 | + # np.column_stack((x, E_interp)) |
| 100 | + # ) |
| 101 | + |
| 102 | + |
| 103 | +# Mesh ------------------------------------------------------------------------ |
| 104 | + |
| 105 | +# Define an elongated cylinder along the y direction |
| 106 | +def cylinder(r, centre, radius): |
| 107 | + if (r[0] - centre[0]) ** 2. + (r[1] - centre[1]) ** 2 <= radius ** 2.: |
| 108 | + return Ms |
| 109 | + else: |
| 110 | + return 0 |
| 111 | + |
| 112 | +# Finite differences mesh |
| 113 | +mesh = CuboidMesh(nx=6, ny=6, nz=35, |
| 114 | + dx=2, dy=2, dz=2, |
| 115 | + unit_length=1e-9) |
| 116 | + |
| 117 | +centre = (np.max(mesh.coordinates[:, 0]) * 0.5, |
| 118 | + np.max(mesh.coordinates[:, 1]) * 0.5) |
| 119 | + |
| 120 | + |
| 121 | +# Prepare simulation ---------------------------------------------------------- |
| 122 | + |
| 123 | +def elongated_part_sim(): |
| 124 | + sim = Sim(mesh) |
| 125 | + sim.Ms = lambda r: cylinder(r, centre, 8) |
| 126 | + sim.add(UniformExchange(A=A)) |
| 127 | + sim.add(UniaxialAnisotropy(Kx, axis=(0, 0, 1))) # Anisotropy along y |
| 128 | + sim.add(Demag()) |
| 129 | + |
| 130 | + return sim |
| 131 | + |
| 132 | + |
| 133 | +# ----------------------------------------------------------------------------- |
| 134 | + |
| 135 | +barriers = {} |
| 136 | + |
| 137 | + |
| 138 | +def test_energy_barrier_cylinder(): |
| 139 | + init_im = [(0, 0, -1), (0, 0.9, 0.1), (0, 0, 1)] |
| 140 | + interp = [8, 8] |
| 141 | + |
| 142 | + for method in ['Geodesic', 'Cartesian']: |
| 143 | + relax_neb(elongated_part_sim(), |
| 144 | + 1e4, 2000, |
| 145 | + 'neb_cylinder_z-axis_{}'.format(method), |
| 146 | + init_im, |
| 147 | + interp, |
| 148 | + save_every=5000, |
| 149 | + method=method |
| 150 | + ) |
| 151 | + |
| 152 | + # Get the energies from the last state |
| 153 | + data = np.loadtxt('neb_cylinder_z-axis_{}_energy.ndt'.format(method))[-1][1:] |
| 154 | + ebarrier = np.abs(np.max(data) - np.min(data)) / (1.602e-19) |
| 155 | + barriers[method] = ebarrier |
| 156 | + |
| 157 | + # assert ebarrier < 0.017 |
| 158 | + # assert ebarrier > 0.005 |
| 159 | + |
| 160 | + for key in barriers.keys(): |
| 161 | + print(key, barriers) |
| 162 | + |
| 163 | + |
| 164 | +if __name__ == '__main__': |
| 165 | + test_energy_barrier_cylinder() |
0 commit comments