Various models of molecular evolution have been proposed that explain how the accumulation of individually neutral mutations can change the functional impact of subsequent mutations. Over time, this phenomenon is expected to lead to fitness drift where, under constant negative selection, the relative fitness of possible sequence states at a position diverges while the fitness of the entire molecule stays the same. Because this is expected to substantially complicate the interpretation of human genomic variants using comparative data, we decided to explore the extent of fitness drift in a model of the evolutionary dynamics of protein structure. This question has recently been addressed in several high profile studies that conflicted in their conclusions. We believe that these studies probably underestimated how much fitness drift occurs by only allowing mutations to go to fixation one at a time, which unrealistically limits the exploration of sequence space by making compensatory substitutions nearly impossible. We have developed a simulation pipeline to overcome this problem and use it to exhaustively explore drift in the ‘epistasis landscape’, which we define as a summary of the fitness impact of all possible pairs of substitutions at all possible pairs of positions. Although this is extremely computationally difficult to evaluate under the molecular mechanics force-field that we employ, it provides a unique window into how the architecture of protein structure may change over evolutionary time that would be nearly impossible to explore experimentally.