Implementing integrator

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 <F | F>
+                force_norm_2 = np.einsum('i,i', force, force)
+                factor = velocity_proj / force_norm_2
+                # Set updated velocity: v = v * (<v|F> / |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)

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

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