We present an approach for computing long-range van der Waals (vdW) interactions between complex molecular systems and arbitrarily shaped macroscopic bodies, melding atomistic treatments of electronic fluctuations based on density functional theory in the former with continuum descriptions of strongly shapedependent electromagnetic fields in the latter, thus capturing many-body and multiple scattering effects to all orders. Such a theory is especially important when considering vdW interactions at mesoscopic scales, i.e., between molecules and structured surfaces with features on the scale of molecular sizes, in which case the finite sizes, complex shapes, and resulting nonlocal electronic excitations of molecules are strongly influenced by electromagnetic retardation and wave effects that depend crucially on the shapes of surrounding macroscopic bodies. We show that these effects together can modify vdW interaction energies and forces, as well as molecular shapes deformed by vdW interactions, by orders of magnitude compared to previous treatments based on Casimir-Polder, nonretarded, or pairwise approximations, which are valid only at macroscopically large or atomic-scale separations or in dilute insulating media, respectively. DOI: 10.1103/PhysRevLett.118.266802 Van der Waals (vdW) interactions play an essential role in noncovalent phenomena throughout biology, chemistry, and condensed-matter physics [1][2][3]. It has long been known that vdW interactions among a system of polarizable atoms are not pairwise additive but instead strongly depend on geometric and material properties [2,4,5]. However, only recently developed theoretical methods have made it possible to account for short-range quantum interactions in addition to long-range many-body screening in molecular ensembles [3,[6][7][8][9][10][11][12][13][14][15], demonstrating that nonlocal many-body effects cannot be captured by simple, pairwise-additive descriptions; these calculations typically neglect electromagnetic retardation effects in molecular systems. Simultaneously, recent theoretical and experimental work has characterized dipolar Casimir-Polder (CP) interactions between macroscopic metallic or dielectric objects and atoms, molecules, or Bose-Einstein condensates, further extending to nonzero temperatures, dynamical situations, and fluctuations in excited states (as in so-called Rydberg atoms) [16][17][18][19][20][21][22][23][24][25]. Yet, while theoretical treatments have thus far accounted for the full electrodynamic response of macroscopic bodies (including retardation), they often treat molecules as point dipoles of some effective bulk permittivities or as collections of noninteracting atomic dipoles, ignoring finite size and other many-body electromagnetic effects.In this Letter, motivated by the aforementioned theoretical developments [1,[16][17][18][24][25][26][27][28], we describe an approach that seamlessly connects atomistic descriptions of large molecules to continuum descriptions of arbitrary macroscopic bodies, characterizing t...