Positive muons implanted in metals tend to stop at interstitial sites that correspond to the maxima of the Coulomb potential energy for electrons in the material. In turns the Coulomb potential is approximated by the Hartree pseudo-potential obtained from the GPAW calculation. A good guess is therefore given by the maxima of this potential.
In this tutorial we obtain the guess in the case of MnSi. The results can be compared with A. Amato et al. [Amato14], who find a muon site at fractional cell coordinates (0.532,0.532,0.532) by DFT calculations and by the analysis of experiments.
Let’s perform the calculation in ASE, starting from the space group of MnSi, 198, and the known Mn and Si coordinates.
from gpaw import GPAW, PW, MethfesselPaxton from ase.spacegroup import crystal from ase.io import write a = 4.55643 mnsi = crystal(['Mn', 'Si'], [(0.1380, 0.1380, 0.1380), (0.84620, 0.84620, 0.84620)], spacegroup=198, cellpar=[a, a, a, 90, 90, 90]) for atom in mnsi: if atom.symbol == 'Mn': atom.magmom = 0.5 mnsi.calc = GPAW(xc='PBE', kpts=(2, 2, 2), mode=PW(800), occupations=MethfesselPaxton(width=0.005), txt='mnsi.txt') mnsi.get_potential_energy() mnsi.calc.write('mnsi.gpw') v = mnsi.calc.get_electrostatic_potential() write('mnsi.cube', mnsi, data=v)
The ASE code outputs a Gaussian cube file, mnsi.cube, with volumetric data of the potential (in eV) that can be visualized.
Getting the maximum¶
One way of identifying the maximum is by the use of an isosurface (or 3d contour surface) at a slightly lower value than the maximum. This can be done by means of an external visualization program, like eg. majavi (see also Plotting iso-surfaces with Mayavi):
$ python3 -m ase.visualize.mlab -c 11.1,13.3 mnsi.cube
The parameters after -c are the potential values for two countour surfaces (the maximum is 13.4 eV).
This allows also secondary (local) minima to be identified.
A simplified procedure to identify the global maximum is the following
# Creates: pot_contour.png from gpaw import restart import matplotlib.pyplot as plt import numpy as np mnsi, calc = restart('mnsi.gpw', txt=None) v = calc.get_electrostatic_potential() a = mnsi.cell[0, 0] n = v.shape x = y = np.linspace(0, a, n, endpoint=False) f = plt.figure() ax = f.add_subplot(111) cax = ax.contour(x, y, v[:, :, n // 2], 100) cbar = f.colorbar(cax) ax.set_xlabel('x (Angstrom)') ax.set_ylabel('y (Angstrom)') ax.set_title('Pseudo-electrostatic Potential') f.savefig('pot_contour.png')
The figure below shows the contour plot of the pseudo-potential in the plane z=2.28 Angstrom containing the maximum
The absolute maximum is at the center of the plot, at (2.28,2.28,2.28), in Angstrom. A local maximum is also visible around (0.6,1.75,2.28), in Angstrom.
In comparing with [Amato14] keep in mind that the present examples has a very reduced number of k points and a low plane wave cutoff energy, just enough to show the right extrema in the shortest CPU time.