Bookkeeping

Here we describe how to implement Algorithm 2 in Fast electrostatic solvers for kinetic Monte Carlo simulations with a ParticleLoop and a PairLoop. The ParticleLoop and PairLoop are the fundamental looping operations of the performance-portable framework for Molecular Dynamics (PPMD) https://github.com/ppmd/ppmd. We use the following definitions for ParticleLoop and PairLoop.

ParticleLoop = loop.ParticleLoopOMP
PairLoop = pairloop.CellByCellOMP

We assume a state A is initialised as follows, as described in propose_with_dats.

# Create a state "A"

A = state.State()

# Set the number of particles
A.npart = N

# Set the domain (must be cubic)
A.domain = domain.BaseDomainHalo(extent=(E,E,E))
# set the boundary condition for the domain (note this does not
# set the domain for the FMMKMC)
A.domain.boundary_condition = domain.BoundaryTypePeriodic()

# Create a PositionDat for particle positions and a ParticleDat
# for particle charge values
A.P = PositionDat(ncomp=3)
A.Q = ParticleDat(ncomp=1)

# Create dats for other desired properties such as global id
A.GID = ParticleDat(ncomp=1, dtype=INT64)

# ParticleDat for ``current_sites``
A.sites = ParticleDat(ncomp=1, dtype=INT64)
# ParticleDat for ``prop_positions``
A.prop_positions = ParticleDat(ncomp=M*3)
# ParticleDat for ``prop_masks``
A.prop_masks = ParticleDat(ncomp=M, dtype=INT64)
# ParticleDat for ``prop_energy_diffs``
A.prop_diffs = ParticleDat(ncomp=M)

# ScalarArray that holds the number of moves per site type
# this example uses a cubic lattice with one site type,
site_max_counts = ScalarArray(ncomp=1, dtype=INT64)

The first step of the algorithm is to loop over all particles and populate the ParticleDat prop_positions with a ParticleLoop. In this example the particles move on a cubic lattice. The set of proposed moves are computed by taking the current particle position and adding an offset for each proposed position. In this example there are M proposed positions for each site (and M corresponding offsets).

# make proposed positions kernel
prop_pos_kernel_src = r'''
// current pos
const double p0 = P.i[0];
const double p1 = P.i[1];
const double p2 = P.i[2];

// reset mask
for (int mx=0 ; mx<M ; mx++){
    MASK.i[mx] = 1;
}

// form proposed positions from offsets
for (int mx=0 ; mx< M ; mx++){
    double n0 = p0 + OA[mx*3 + 0];
    double n1 = p1 + OA[mx*3 + 1];
    double n2 = p2 + OA[mx*3 + 2];

    // wrap positions into domain
    if ( n0 < LOWER ) { n0 += EXTENT; }
    if ( n1 < LOWER ) { n1 += EXTENT; }
    if ( n2 < LOWER ) { n2 += EXTENT; }
    if ( n0 > UPPER ) { n0 -= EXTENT; }
    if ( n1 > UPPER ) { n1 -= EXTENT; }
    if ( n2 > UPPER ) { n2 -= EXTENT; }

    PP.i[mx*3 + 0] = n0;
    PP.i[mx*3 + 1] = n1;
    PP.i[mx*3 + 2] = n2;
}
'''

# create a kernel from the kernel source
prop_pos_kernel = kernel.Kernel(
    'prop_pos_kernel',
    prop_pos_kernel_src,
    constants=(
        Constant('M', M),
        Constant('LOWER', -0.5 * E),
        Constant('UPPER', 0.5 * E),
        Constant('EXTENT', E)
    )
)

# create a ParticleLoop from the kernel
prop_pos = ParticleLoop(
    kernel=prop_pos_kernel,
    dat_dict={
        'P'     : A.P(READ),
        'PP'    : A.prop_positions(WRITE),
        'OA'    : offsets_sa(READ),
        'MASK'  : A.prop_masks(WRITE)
    }
)

The ParticleLoop prop_pos is executed by calling

prop_pos.execute()

Now that the proposed positions are stored in prop_pos we remove proposed moves that would result in particle overlap by masking off “bad” moves. We discover potential overlaps with the PairLoop operation. For each particle we loop over neighbouring particles within max_move and check for an overlap with the proposed moves. If an overlap is detected then the corresponding move is removed by setting prop_masks to a value less than 1 for that move. We allow for any value less than 1 such that recombination events can be identified, the current FMM-KMC implementation does not support recombination.

We provide now the example kernel to mask off proposed moves that would cause an overlap. It is important to note that PPMD does not wrap particle positions around the periodic boundary. To be explicit: if a particle has a \(x\) position of \(-0.5E + \delta\) then it will be exposed with a \(x\) position of \(0.5E + \delta\) when observed over the periodic boundary in a PairLoop. This convention avoids the use of conditionals in PairLoops for inter-particle interactions, which is the primary use case of the PPMD framework.

Bearing this in mind, we recompute the proposed moves for each particle for each neighbour, without wrapping around the boundary, and then compare particle positions.

# make exclude kernel

exclude_kernel_src = r'''
// current position of particle i
const double p0 = P.i[0];
const double p1 = P.i[1];
const double p2 = P.i[2];

// read particle j position
const double pj0 = P.j[0];
const double pj1 = P.j[1];
const double pj2 = P.j[2];

// check each proposed position
for (int mx=0 ; mx< M ; mx++){

    // form the proposed position without wrapping
    double n0 = p0 + OA[mx*3 + 0];
    double n1 = p1 + OA[mx*3 + 1];
    double n2 = p2 + OA[mx*3 + 2];

    // difference to particle j position
    const double d0 = pj0 - n0;
    const double d1 = pj1 - n1;
    const double d2 = pj2 - n2;

    // if they overlap, mask off the position
    const double r2 = d0*d0 + d1*d1 + d2*d2;

    // using a ternary operator aids autovectorisation.
    MASK.i[mx] = (r2 < TOL) ? 0 : MASK.i[mx];
}
'''

# create kernel from kernel source
exclude_kernel = kernel.Kernel(
    'exclude_kernel',
    exclude_kernel_src,
    constants=(
        Constant('M', M),
        Constant('TOL', 0.01)
    )
)

# create pairloop from kernel
exclude = PairLoop(
    kernel=exclude_kernel,
    dat_dict={
        'P'     : A.P(READ),
        'OA'    : offsets_sa(READ),
        'MASK'  : A.prop_masks(WRITE)
    },
    shell_cutoff = max_move
)

As for the ParticleLoop, the PairLoop exclude can be executed with

exclude.execute()