diff --git a/fidimag/common/chain_method_integrators.py b/fidimag/common/chain_method_integrators.py index 01c379bd..1aa6f10a 100644 --- a/fidimag/common/chain_method_integrators.py +++ b/fidimag/common/chain_method_integrators.py @@ -10,8 +10,8 @@ class StepIntegrator(BaseIntegrator): """ - A simple integrator where spins are normalised at every inetegrator step - Integrator options are Euler and RK4 + A simple integrator where spins are normalised at every integrator step. + Integrator options are `euler` and `rk4` """ def __init__(self, spins, rhs_fun, step="euler", stepsize=1e-15): super(StepIntegrator, self).__init__(spins, rhs_fun) @@ -36,13 +36,13 @@ def set_options(self, rtol=1e-8, atol=1e-8): def set_step(self, step): step_choices = {'euler': euler_step, 'rk4': runge_kutta_step} if step not in step_choices: - raise NotImplemented("step must be euler or rk4") + raise NotImplementedError("step must be euler or rk4") self._step = step_choices[step] class VerletIntegrator(BaseIntegrator): - """ - A quick Verlet integration in Cartesian coordinates + """A quick Verlet integration in Cartesian coordinates + See: J. Chem. Theory Comput., 2017, 13 (7), pp 3250–3259 """ def __init__(self, band, forces, rhs_fun, n_images, n_dofs_image, @@ -142,6 +142,98 @@ def _step(self, t, y, h, f): return tp +# ----------------------------------------------------------------------------- +# NEBM integrator for F-S algorithm + + +class FSIntegrator(BaseIntegrator): + """A step integrator considering the action of the band + """ + def __init__(self, band, forces, action, rhs_fun, n_images, n_dofs_image, + max_steps=1000): + super(FSIntegrator, self).__init__(band, rhs_fun) + + self.i_step = 0 + self.n_images = n_images + self.n_dofs_image = n_dofs_image + self.forces_prev = np.zeros_like(band).reshape(n_images, -1) + # self.G : + self.forces = forces + self.max_steps = max_steps + + def run_until(self, t): + pass + + def run_for(self, n_steps): + # while abs(self.i_step - steps) > EPSILON: + st = 1 + while st < n_steps: + self.i_step = self._step(self.t, self.y, self.stepsize, self.rhs) + + self.rhs_evals_nb += 0 + + if self.i_step > self.max_steps: + break + + st += 1 + return 0 + + def set_options(self): + pass + + def _step(self, t, y, h, f): + """ + """ + + f(t, y) + force_images = self.forces + force_images.shape = (self.n_images, -1) + # In this case f represents the force: a = dy/dt = f/m + # * self.m_inv[:, np.newaxis] + y.shape = (self.n_images, -1) + # print(force_images[2]) + + # Loop through every image in the band + for i in range(1, self.n_images - 1): + + force = force_images[i] + velocity = self.velocity[i] + velocity_new = self.velocity_new[i] + + # Update coordinates from Newton eq, which uses the "acceleration" + # At step 0 velocity is zero + y[i][:] = y[i] + h * (velocity + (h / (2 * self.mass)) * force) + + # Update the velocity from a mean with the prev step + # (Velocity Verlet) + velocity[:] = velocity_new[:] + (h / (2 * self.mass)) * (self.forces_prev[i] + force) + + # Project the force of the image into the velocity vector: < v | F > + velocity_proj = np.einsum('i,i', force, velocity) + + # Set velocity to zero when the proj = 0 + if velocity_proj <= 0: + velocity_new[:] = 0.0 + else: + # Norm of the force squared + force_norm_2 = np.einsum('i,i', force, force) + factor = velocity_proj / force_norm_2 + # Set updated velocity: v = v * ( / |F|^2) + velocity_new[:] = factor * force + + # New velocity from Newton equation (old Verlet) + # velocity[:] = velocity_new + (h / self.mass) * force + + # Store the force for the Velocity Verlet algorithm + self.forces_prev[:] = force_images + + y.shape = (-1,) + force_images.shape = (-1,) + normalise_spins(y) + tp = t + h + return tp + + def normalise_spins(y): # Normalise an array of spins y with 3 * N elements y.shape = (-1, 3) diff --git a/fidimag/common/nebm_FS.py b/fidimag/common/nebm_FS.py index 04bdbabf..20fa907b 100644 --- a/fidimag/common/nebm_FS.py +++ b/fidimag/common/nebm_FS.py @@ -148,6 +148,21 @@ def initialise_energies(self): self.energies[i] = self.sim.compute_energy() self.band = self.band.reshape(-1) + def initialise_integrator(self): + self.t = 0 + self.iterations = 0 + self.ode_count = 1 + + self.integrator = FSIntegrator(self.band, # y + self.G, # forces + self.step_RHS, + self.n_images, + self.n_dofs_image, + stepsize=1e-4) + self.integrator.set_options() + # In Verlet algorithm we only use the total force G and not YxYxG: + self._llg_evolve = False + def generate_initial_band(self, method='linear'): """ method :: linear, rotation diff --git a/fidimag/common/string_method.py b/fidimag/common/string_method.py index d3e74cbf..f7520d18 100644 --- a/fidimag/common/string_method.py +++ b/fidimag/common/string_method.py @@ -252,9 +252,7 @@ def compute_distances(self): self.path_distances[1:] /= self.path_distances[-1] def step_RHS(self, t, y): - """ - Use Step integrators from the chain_method_integrators library - + """Use Step integrators from the `chain_method_integrators` library """ self.ode_count += 1 @@ -282,9 +280,9 @@ def Sundials_RHS(self, t, y, ydot): we are solving dy/dt = 0 WARNING: The variable step integrator from Sundials does not work well - with the StringMethod, making the algorithm overshoot the solutions for - large time steps and driving the images toward the extrema images. We - could poossibly fix this by redefining the positions of the images + with the `StringMethod`, making the algorithm overshoot the solutions + for large time steps and driving the images toward the extrema images. + We could possibly fix this by redefining the positions of the images after every integrator step, instead of redefining them after a certain number of steps, however this requires to tune the Sundials Python wrapper. In addition, we would need to check the stopping criteria of