int change_volume_and_shape() {
    // Calculate current volume
    double volume = box[0] * box[1] * box[2];
    double scale_factor = 1.;

    //Select a random axis
    int random_axis = rnd.RandomIntAB(0, NDIM - 1);
    double old_axis_length = box[random_axis];
    box[random_axis] += 3 * rnd.RandomMinusOneOne();
    double new_volume = box[0] * box[1] * box[2];

    //Acceptance criterion for volume change move
    double prob = exp(-beta * pressure * (new_volume - volume) + (n_particles + 1.) * log(new_volume / volume));
    if (rnd.RandomZeroOne() > prob) {
        box[random_axis] = old_axis_length;
        return 0;
    }

    //Scale the particle coordinates and the box
    scale_factor = box[random_axis] / old_axis_length;
    for (int n = 0; n < n_particles; ++n) {
        r[n][random_axis] *= scale_factor;
    }
    if (isThereAnyOverlap()) {
        //Reject move and reset to original packing fraction and return
        //Maybe copy old values instead? The multiplied values differ ~10^-14...
        scale_factor = old_axis_length / box[random_axis];
        for (int n = 0; n < n_particles; ++n) {
            r[n][random_axis] *= scale_factor;
        }
        box[random_axis] = old_axis_length;

        return 0;
    }
    return 1;
}