Next: Numerical Results Up: Determining Pair Interactions from Previous: The Model

Numerical Implementation

In this section we describe our numerical implementation. Simulations were carried out on a square region with 7250 randomly placed pins. Having constructed a random pinscape, we placed N = 65 flux lines at random and solved numerically for stable equilibrium configurations of Eq. (5) using a three stage procedure: (i) Starting from a random initial configuration, the vortex system is allowed to relax dynamically according to

 (9)

until the magnitudes of all forces on the RHS are less than a prescribed target value. (ii) Once this has been achieved, the resulting configuration is used as an initial guess in a Newton-Raphson algorithm to find a solution to Eq. (5). Solutions found this way are not guaranteed to be stable with respect to small perturbations of the particle locations. Therefore, (iii), a linear stability analysis is carried out to check for stability. If the configuration obtained this way does not turn out to be stable, stage (i) is resumed with a reduced target value. The whole procedure is repeated until a stable configuration has been found. The combination of a molecular dynamics algorithm with a Newton-Raphson scheme, as described above, turns out to be far more computationally efficient than using molecular dynamics alone [9].

Lengths were chosen so that the square corresponds to a m2 field of view in a magnetic field of roughly 50 Gauss. With this choice, all lengths will be reported in m. The London penetration depth was taken to be m. These choices resemble the experimental conditions in Ref. [3]. The parameters varied in the simulations were the range of the pinning forces rp, with values m, and the distribution of pinning strengths , cf. Eq. (4). We considered the following distributions:

 (10)

for identical pinning strengths, and

 (11)

for stretched-exponentially distributed pinning strengths, with the constants c1 and c2 given by , and , where is the Gamma function. In particular, the distributions used were (half-gaussian), (exponential), and (stretched exponential). The mean of the above distributions is given by Vo. The unit of energy was chosen such that

 (12)

This choice is consistent with experimental observations [3], and ensures that the maximum force exerted by a pin is independent of the range of the pinning potential.

The numerical integration in stage (i) was carried out using a 4th order Runge-Kutta method [10] with adaptive step size and an initial force resolution of 10-6 in the above units with subsequent decrements by factors of 10, as needed. Algorithms of the software libraries LINPACK and EISPACK [11] were used in putting together stages (ii) and (iii). The maximum residual force on any particle in a stable configuration was found to be smaller than .

Free boundary conditions were employed, and flux lines leaving the square (leaks) were put back in random locations. Typically, only a few leaks occurred in a given run.

In order to speed up computations further, several additional measures were taken: If a stable configuration was not found after a fixed number of integration steps (5000), a new initial configuration was generated. Initial configurations were such that flux lines were randomly positioned, but no two were allowed to be closer initially than m. These measures were found to reduce the computation time without changing the results significantly. The particle-pin interaction, Eq. (4), was cut off at a distance of 4rp. Prior to each run, the area was subdivided into a mesh of squares and a look-up table was generated providing each square with the location and strength of pins within the effective interaction distance 4rp. This eliminates the need to search all pinning sites for each computation of the pinning forces.

Next: Numerical Results Up: Determining Pair Interactions from Previous: The Model
David G. Grier
1998-11-19