This tutorial is suited to those interested in carrying out DFT+U calculations on a system with multiple transitionmetal sites. If, instead, you are only concerned with singlesite complexes, check out the original tutorial on calculating the Hubbard U. You may wish to revisit the singlesite U tutorial for more background and instruction.
Background: The DFT+U approach in QuantumESPRESSO is equipped to not only treat singlesite complexes but multiple sites as well. Typically, the extension to multiple sites is important where there is more than one unique transitionmetal in the system or, in rare cases, where an onsite U on ligand atoms is useful in addition to transitionmetals.
The multiplesite, linearresponse U is determined from a matrix of linearresponse functions of the form[1]:
where the relationship between the perturbation of electrons on site J and the response on site I yields the response function denoted by an index IJ. In order to obtain U for multiple sites, we now invert the converged and bare matrices[1]. The value of U for site I is given by:
where the offdiagonal elements indirectly influence the value of onsite Uthrough the matrix inversion.
Extensions to this approach in which we treat multiple manifolds on the
same site  e.g. both 3d and 4s electrons on a transition metal[2]  or directly treat intersite coupling from offdiagonal terms[3] are both possible but neither are currently available in a public version of QuantumESPRESSO. These more advanced topics may be discussed in the future.
Instruction: As was the case for singlesite, linearresponse U calculation, the extension to a matrix of response functions adds little additional computational cost. We extend the previous case of ground state sextet MnO to include a calculation of linearresponse Hubbard U values for both the Mn 3d manifold and the O 2p manifold.
We determine the linearresponse U for a multiplesite complex as follows:

1.Obtain single point energy at zero alpha and store wavefunctions.

2.Starting from 1, obtain new single point energy at several values ofnonzero alpha on site J (e.g. 0.08 to 0.08) with tight convergence criteria.

3.Collect occupations of all sites from first iteration of 2 for bare response and last iteration of 2 for converged response.

4.Calculate response functions from linear regression of all relationships obtained in 3 via linear regression.

5.Repeat 24 by applying alpha sequentially to each additional nonJ site.

6.Invert the matrices and subtract to obtain each U from diagonal elements.
After jobrun.py carries out calculations, the linregress.py script calculates response functions and inverts the matrices. If matplotlib is enabled, we can visualize our results on MnO to obtain the following:
The recalculated U on Mn 3d is now 2.9 eV, reduced from the value of about 3.5 eV obtained when only calculating U on a single site. These two treatments produce slightly different values of U because the multiplesite case includes the effects of O 2p.
Note that we obtain the inverse in 6 routinely for the 2x2 case. For larger matrices, the script uses the scipy module linalg. If linregress.py detects that scipy is not installed, it will write the linearresponse matrices and you will need to invert the matrices with the software package of your choosing.
Summary: The tutorial files, provided also as a zipped archive here, are:

•jobrun.py — skeleton script generates input files and runs jobs.

•atreader.py — parses your xyz file to generate run parameters.

•variables.py — you should change these job and cluster variables!

•libraries.py — dictionary of parameters that should not be changed.

•linregress.py — calculates matrix of U values from jobrun.py results.

•other files — pseudopotentials, coordinates for MnO, a readme file.
Advanced note: If you try to run DFT+U on an unconventional element, the code may exit with the error “Pseudopotential not yet inserted”. In order to circumvent this issue, you will want to modify the files set_hubbard_l.f90 and tabd.f90 in the PW folder of your QuantumESPRESSO root directory.
I hope that this tutorial has helped you to better understand how to calculate the linearresponse and selfconsistent Hubbard U for multiplesite transition metal complexes. Please email me if you have any additional questions not answered here!
References:
[1] M. Cococcioni and S. de Gironcoli. Physical Review B 71, 035105 (2005) and references therein for more details.
[2] H. J. Kulik, M. Cococcioni, D. A. Scherlis, and N. Marzari. Physical Review Letters 97, 103001 (2006).
[3] H. J. Kulik and N. Marzari. Journal of Chemical Physics 134, 094103 (2011).