We present a first-principles scheme for computing the magnetoelectric response of multiferroics. We apply our method to BiFeO3 (BFO) and related compounds in which Fe is substituted by other magnetic species. We show that under certain relevant conditions -i.e., in absence of incommensurate spin modulation, as in BFO thin films and some BFO-based solid solutions -these materials display a large linear magnetoelectric response. Our calculations reveal the atomistic origin of the coupling and allow us to identify the most promising strategies to enhance it.PACS numbers: 75.80.+q, 71.15.Mb Magnetoelectric (ME) multiferroics present coupled electric and magnetic orders [1], which may allow the development of a variety of magnetic devices, from memories to spin filters, whose behavior would be switchable by application of a voltage. The experimental search for robust room-temperature (T room ) multiferroics is proving a major challenge. The first-principles contribution to this effort is quickly increasing, as shown by recent predictions of new materials [2,3,4] and novel ME coupling mechanisms [5,6,7]. Yet, many key issues remain to be addressed theoretically, as even the ME response of the most promising systems (e.g., thin films of BiFeO 3 [8]) still needs to be characterized and understood in detail. Here we describe a method for ab initio computations of the ME response of multiferroics, tackling the strainmediated contributions that typically exist in these materials. We apply our method to BiFeO 3 and related compounds.Methodology.-The response of a linear magnetoelectric that is not multiferroic -and is thus paraelectriccan be split in two contributions [9]: a purely electronic one and an ionic one. The ionic part accounts for the structural response to an electric (magnetic) field and the resulting change of magnetization (polarization). When the linear magnetoelectric is multiferroic -and thus ferroelectric -there is a third contribution associated to the piezoelectric response to the electric field. Indeed, it can be shown that all commensurate [11] multiferroics must be piezomagnetic, as the symmetries that preclude piezomagnetism are broken in such systems [12]. Hence, the combination of piezoelectricity and piezomagnetism in multiferroics results in a strain-mediated contribution to the linear ME response.In order to model such effects, we have generalized to the magnetoelectric case the formalism that Wu, Vanderbilt, and Hamann [13] (WVH in the following) introduced for a systematic ab initio treatment of dielectric and piezoelectric responses. Let us consider how the energy of a multiferroic crystal varies as a function of the following perturbations: displacements u m of atoms away from their equilibrium positions, homogeneous strains η j , and applied electric (E α ) and magnetic (H µ ) fields. Here, m is a composite label that runs over the atoms in the unit cell and displacement directions, j runs from 1 to 6 in Voigt notation, and α and µ are Cartesian directions.We write the energy per unde...