in DFT program such as VASP, this can be simply calculated.
For example, set all spin directions along 001 and 111, then compare two total energies. Of course, you need apply noncollinear spin polarization and SOC as well.
It is my understanding that you cannot do that with VASP 5.2.1. I don't know about newer versions. I believe you can only apply initial magnetic moments (a guess) and/or external magnetic field (for NMR). Perhaps you can loop with varying external field (may be worth looking into) for a hysteresis plot, but I wouldn't count on it. Maybe it is also worth looking at wannier90 (a post processing software), but no promises. I'm just throwing some ideas out there, but I think what you want cannot be done currently with VASP.
It's not that simple. The magnetic anisotropy energy is just one ingredient. Another one is magnetostriction (the strain-dependence of magnetic anisotropy) and the exchange constant and finally the real structure of your system.
Initially, you have magnetic domains separated by Bloch walls. The domain structure occurs in order to minimize the magnetic polarization field energy (having macroscopic multipoles that vanish or are as small as possible). External fields change the domain structure. Magnetic moments can be rotated out of "easy axis" directions. Domains with favorably aligned magnetic moments grow, the others shrink (and finally one macroscopic single domain remains). Going back to zero external field the reality is that despite the tendency to go back to the initial minimum-energy domain structure there remains a finite magnetization. The reason is pinning of Bloch walls at defects (e.g. grain boundaries, dislocations, point defects) which prevents this simple back-relaxation. Just some extra "force", e.g. due to an external magnetic field, can start to move out some of the domain walls from their pinned positions and at a certain (coercive) field we reach again a vanishing macroscopic magnetization (and finally it goes in the opposite direction). By the way: A strong mechanical impact or heating can also unpin the Bloch walls and demagnetize the system.
All this finally defines your hysteresis curve and quantities like remanence and coercivity. The crucial point is: VASP could help you to determine material parameters like the magnetic anisotropy energy (numerically critical), maybe even magnetostriction constants (numerically extremely critical) or magnetic exchange energy. But it cannot provide the remaining parameters. You have to put that then into a macroscopic continuum theory ("micromagnetism"). There you need to know the shape and structure of your magnet, what kind of defects and how many of them it possesses and which strain fields are related to these defects (in order to be able to determine Bloch wall pinning forces). So, it depends on the individual sample and results may vary by orders of magnitude from sample to sample ...
There is only one really simple case: Very small so-called single-domain particles. There you do not need to bother about Bloch wall pinning etc. (because there are no domain walls), there it is pure rotation of the magnetic moments against magnetic anisotropy. Such single-domain particles (many of them embedded in a usually non-magnetic matrix and pre-oriented through a very strong magnetic field during the fabrication process) are used for real "permanent" magnets since even serious mechanical impacts cannot demagnetize them. If you had something like that in mind then the simple reduction to magnetic anisotropy coercivity works. In any other case it's impossible to rely on a simple ab-initio study and there will be no universal answer (it's a wide range of possible results depending on the defect densities and defect distribution and it requires a continuum theory simulation on top -- you could just calculate the material parameters entering the continuum theory equations ...).
Maybe you could find more details and insights about all that in the book of Helmut Kronmueller and Manfred Faehnle titled "Micromagnetism and the Microstructure of Ferromagnetic Solids" by Cambridge Univ. Press (2003).
You are right, my question was from what you call the single-domain particles point of view. Of course for a realistic material micromagnetism must be used!
However, I believe that the magnetic anisotropy is related to the CHANGE of the coercive field, more than with its absolute value.
I was expecting to obtain from DFT calculation, something like the coercive field for a domain, which could be used as an input for a micromagnetic algorithm.
I think that a realistic material is the "sum" (better a "conjunction") of domains, where for each one, i can get accurate information from DFT algorithm.
Also I hope to obtain with DFT something like an upper bond for the coercive field.
Not directly for the coercive field. You can calculate the anisotropy energy. For a single-domain particle you can then do a very simplified micromagnetism (you have the interaction energy between magnetization and external field and the anisotropy energy (and if there is a surface also the magnetic polarization field energy); you can then look for the minimum-energy configuration, i.e., equilibrium direction of the magnetic moments and where they "switch" into/towards the magnetic field direction). From that I would guess that anisotropy and coercive field should be proportional (but may depend on the orientation of the magnetic field relative to the "easy axis"; "off-axis" I would expect that the magnetization approaches the magnetic field direction asymptotically only, just "on-axis" it is a hard , e.g. 180 degree, switch at a certain critical field -- at T=0K). Of course, take into account that finite temperature effects ("thermal excitation") lead to smaller coercive fields (also by the temperature dependence of the anisotropy energy). So the upper bound should be the T=0K result ... .
However, this really only works for a true single-domain particle. If several domains are involved then the motion of domain walls (growth and shrinkage of domains) is the important process (formally no "barrier" involved, in reality small barriers due to domain wall pinning at defects -- this requires definitely also the inclusion of the magnetic polarization field energy, i.e. of surfaces/crystal shape), not a "switching" of the magnetization direction (reduction of an anisotropy barrier by an external field). The latter might somehow still define an upper bound but I think an unrealistic one (I guess the true coercive field based on domain wall motion should be much smaller ...). But what I understand is that you just want to skip all these domain wall motion processes(?) and I think it will not give any reliable or useful result if you do so. The single-domain limit is a bit too simple I think and it gives an unrealistic (much too large) upper bound(?).
Our steel has a coercive field near 700 A/m which is near 3/2 saturation magnetostriction times residual stress divided by saturation induction, the pinning field. We find 700 A/m to be about 50 A/m less than the pinning field possibly due to the viscous field in the one second range. Many authors use the 100 axis magnetostriction present within each domain to calculate the coercive field, and this is 50 per cent higher for our steel.