Houdini’s Grain solver operates on Position-Based Dynamics (PBD) principles, which decouple collision resolution from velocity integration. Instead of calculating forces and accelerating particles, the solver directly modifies particle coordinates to satisfy spatial constraints, then reconstructs velocities from the resulting positional delta.

The foundational step involves detecting overlapping particles and pushing them apart. Assuming a uniform particle radius, the correction vector accumulates based on penetration depth:

float search_rad = @pscale * 2.0;
int neighbors[] = nearpoints(0, @P, search_rad, 50);

vector total_shift = {0, 0, 0};
int valid_hits = 0;

foreach (int idx; neighbors) {
    if (idx == @ptnum) continue;
    
    vector neighbor_pos = point(0, "P", idx);
    float dist = distance(@P, neighbor_pos);
    
    if (dist < search_rad && dist > 1e-5) {
        vector dir = normalize(@P - neighbor_pos);
        float penetration = (search_rad - dist) * 0.5;
        total_shift += dir * penetration;
        valid_hits++;
    }
}

if (valid_hits > 0) {
    @P += total_shift / valid_hits;
}

When particle mass varies, the displacement must be distributed proportionally. Heavier particles move less, while lighter ones absorb more of the correction:

float search_rad = @pscale * 2.0;
int neighbors[] = nearpoints(0, @P, search_rad, 50);

vector weighted_shift = {0, 0, 0};
int count = 0;

foreach (int idx; neighbors) {
    if (idx == @ptnum) continue;
    
    vector npos = point(0, "P", idx);
    float dist = distance(@P, npos);
    
    if (dist < search_rad && dist > 1e-5) {
        float n_mass = point(0, "mass", idx);
        float influence = n_mass / (@mass + n_mass);
        vector push_dir = normalize(@P - npos);
        float overlap = search_rad - dist;
        
        weighted_shift += push_dir * overlap * influence;
        count++;
    }
}

if (count > 0) {
    @P += weighted_shift / count;
}

Once positions are resolved, velocities are derived by comparing the current frame’s coordinates against the stored previous positions. This bypasses traditional force integration:

if (i@has_pprevious) {
    @v = (@P - v@pprevious) / @TimeInc;
}

The POP Grain node automatically flags particles with an ispbd attribute set to 1. This signal instructs the upstream POP Solver to skip its native position integration step, ensuring that all spatial adjustments occur exclusively within the PBD constraint loop.

A critical convention in Houdini’s grain system involves the mass attribute. Assigning a value of 0 does not represent a massless particle; instead, it signals infinite mass. In the mass-weighted correction formula, as @mass approaches infinity, the influence ratio n_mass / (@mass + n_mass) converges to zero. Consequently, the particle becomes completely immovable during collision resolution, acting as a static collider or anchor.

When implementing this behavior via VEX, attribute initialization requires careful handling due to how Houdini compiles wrangle contexts. Declaring the attribute with a default value across the entire geometry is necesary before applying conditional overrides:

// Mandatory initialization to prevent undefined attribute errors
f@mass = 1.0; 

// Conditional override for specific groups
if (inpointgroup(0, "static_anchors", @ptnum)) {
    @mass = 0.0;
}

Without the initial f@mass = 1.0; declaration, the solver may fail to recognize the attribute scope, leading to inconsistent simulation states. Once applied, particles in the target group remain fixed while surrounding dynamic grains collapse and flow around them under gravity.