Flexible harmonic boundary conditions have been proposed by Sinclair in the 1970s in order to overcome spurious effects on atomistic problems due to fixed boundaries. To date this method has never been applied to problems beyond isolated defects, such as dislocations, because it involves dense boundary matrices which become quickly unsuitable to larger problems due to their vast memory requirements. In order to apply the method for larger systems, e.g. arrangements of defects, we propose an implicit approximate representation using hierarchical matrices which have proven efficiency in the context of boundary integral equations while preserving overall accuracy. Despite its simplicity, Sinclair’s staggered method converges rather slowly if the approximate far-field harmonic response and the true nonlinear atomic response differ considerably. Starting from Sinclair’s iteration equation for the harmonic displacements, we derive a discrete variant of the well-known boundary element method (BEM) for the exterior balance equation which is then combined with the fully atomistic problem. To solve the coupled problem we propose a monolithic Newton–Krylov scheme which iterates simultaneously on all unknowns. We outline the superior performance of this method in comparison to other existing methods and to classical clamped boundary conditions with numerical examples. Further, we present guidelines for an efficient implementation into existing molecular dynamics codes. •Efficient implementation of flexible boundary conditions for atomistic problems.•Discrete boundary element method based on lattice Green functions.•Hierarchical matrix approximation of discrete boundary element matrices.•Efficient monolithic solver for force-based atomistic/continuum coupling methods.•Considerably superior performance of flexible over clamped boundary conditions.