(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
m^{2} 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 *r*_{p}, with values
m,
and the distribution of pinning
strengths ,
*cf.* Eq. (4).
We considered the following distributions:

for

for stretched-exponentially distributed pinning strengths, with the constants

(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 4^{th}
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 4*r*_{p}. 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 4*r*_{p}. This
eliminates the need to search all pinning sites for each computation of
the pinning forces.