tar -zxvf examples_pseudo.tgz cd examples_pseudo/
Example of pseudopotential generation: Silicon
Step 1: All-electron calculationGo to directory "step1":
cd step1/The file "si_ae.in" contains in the &input namelist the data for an all-electron ("iswitch=1") calculation for Si (Z=14, variable "zed"; you may specify "atom='Si'" instead) in the configuration 3s23p2 (this is the ground state), using LDA (Perdew-Zunger), nonrelativistic ("rel=0") calculation: see the &input namelist. Run the ld1.x code with "si_ae.in" as ionput, for instance as follows:
$espresso_dir/bin/ld1.x < si_ae.in > si_ae.outAll-electron Kohn-Sham states are written to file "ld1.wfc" and can be plotted using "gnuplot", or any other plotting program. Column 1 contains the radial grid, columns 2-6 the radial wavefunctions (in reversed order! the last wavefunction, 3p in this case, is the first column). In order to plot 3s,3p states, for instance:
gnuplot> plot [0:5] 'ld1.wfc' u 1:3 w l gnuplot> replot 'ld1.wfc' u 1:2 w lA better plot can be obtained in file "gnuplot.ps" (remove or comment the three top lines to see the plot at the terminal) using the script "plotwfc.gnu":
gnuplot> load 'plotwfc.gnu'The subdirectory "reference/" contains the reference output for comparison. Does the output look as you expect?
Step 2: Calculation of logarithmic derivativesGo to directory "step2":
cd ../step2/We need to specify the radius "rlderiv" at which "nld" logarithmic derivatives are calculated, in an energy range delimited by "eminld" and "emaxld", on a grid with "deld" spacing (energies in Ry). Run the ld1.x code with input data "si_ae_ld.in". Since we set nld=3, the s,p,d logarithmic derivatives are calculated and written to file "ld1.dlog" (file format: energies in the first column, s,p,d,...logarithmic derivatives in the following columns). Plot them with gnuplot:
gnuplot> plot [-11:2] 'ld1.dlog' u 1:4 w l gnuplot> replot 'ld1.dlog' u 1:3 w l gnuplot> replot 'ld1.dlog' u 1:2 w lDon't be scared by the wild oscillations: it is just the state going through 0. For a better plot in file "gnuplot.ps", use the script "plotld.gnu":
gnuplot> load 'plotld.gnu'Script "plotld_sq.gnu" will plot the same, plus a square indicating the energy range of interest for valence states. This is the region that a pseudopotential should be able to reproduce.
Step 3: generation of the pseudopotentialLet us now generate a pseudopotential for Si. Move to directory "step3". File "si_ps.in" contains the data needed to generate a pseudopotential ("iswitch=3"), starting from the same reference configuration as above. In namelist &inputp one has to provide the channel chosen as local potential (the 3d state: "lloc=2"), the type of pseudopotential (one projector per channel, "pseudotype=1"), the number of valence electrons ("zval=4"), the output pseudopotential file. The following lines contain a list of states to be pseudized, with corresponding pseudization energies and matching radii (in this case, rc=2.4 for s,p,d). The local channel must be the last! The Rabe-Rappe-Kaxiras-Joannopoulos pseudization technique is used.
Note the special treatment of the 3d state: since it is not bound on the ground state, it is flagged by a (bogus) negative value of the charge in the all-electron configuration, and it is pseudized at energy -0.10Ry, chosen more or less arbitrarily in the typical energy range for valence electrons; 3s and 3p are instead pseudized at the Kohn-Sham energy of the bound state (this is the meaning of "0.00" as pseudization energy).Run the code with the "si_ps.in" input. The output pseudopotential is written to file Si.rrkj3.UPF.
Note that the suffix ".UPF" implies that the pseudopotential is written in the separable form even if it is generated in the semilocal form.
Check that the one-electron levels are the same for all-electron and pseudopotential calculations. Compare the all-electron and pseudo orbitals: "ld1.wfc" and "ld1ps.wfc", respectively. For a nice plot, use the script "plotwfc.gnu" (arrows point to the matching radii).
Note that the relative sign of all-electron and pseudo-orbitals is arbitrary! this is why the pseudo 3p state is plotted with reversed sign.Do the same for the logarithmic derivatives: "ld1.dlog" and "ld1ps.dlog", respectively. Important: choose an energy range that is is relevant for valence electrons. You can use the script "plotld.gnu"; arrows locate the pseudization energy (where all-electron and pseudo logarithmic derivatives are equal by construction). Nice, isn't it? but this is a simple case, and we are looking at the reference configuration only.
Step 4: test of the pseudopotentialAs a strict minimun, one has to test for transferability on a few atomic configurations that differ from the reference one (the one used to produce the pseudopotential). Move to directory "step4" and copy there the pseudopotential file obtained in the previous step:
cd ../step4/ cp ../step3/Si.rrkj3.UPFFile "si_test.in" contains data for testing ("iswitch=2") over 5 ("nconf") atomic configurations supplied in namelist &test.
Note that fractionary occupation numbers are allowed: mathematically, it is a well defined problem; physically, you may think that you are mimicking an atom in a molecule or a solid.Run the code with input data "si_test.in". A summary of the results is written to file "ld1.test". All-electron and pseudo-energy differences are referred to the first specified configuration, in this case the ground state.
Note that absolute pseudo-energies have no physical meaning! they depend on the specific pseudopotential. Only energy differences are significant.You may want to verify that the all-electron and pseudo-orbitals vary, not by much but by a visible amount, from a configuration to another: see script "diffwfc.gnu".
Step 5: relativistic pseudopotentialFor Si, relativistic effects are small and can be usually neglected. For heavier atoms, however, you may need to take them into account. The dominant relativistic corrections are easily accounted for as perturbations by the so-called scalar-relativistic approximation: just set "rel=1". More accurate calculations require the solution of the Dirac-Kohn-Sham equation. The l>0 orbitals are classified by a quantum number j=l-1/2 and j=l+1/2. As a consequence the pseudopotential for l>0 channels is also split into j=l-1/2 and j=l+1/2 components. Move to directory "step5". File "si_ps_rel.in" contains data for the generation of a relativistic pseudopotential. Note option "rel=2" (fully relativistic) and the presence of the j quantum number in the list of states to be pseudized. Also note that you need both j-components for the state chosen as local (3d in this case). Run the code, examine the output, perform the usual checks.
Step 6: Ultrasoft pseudopotentialUltrasoft (Vanderbilt) pseudopotential allow a drastic reduction of the basis set dimension in plane-wave calculations, at the expense of a more complicated problem to solve. For Si it is not really needed to generate an ultrasoft pseudopotential: notice the "estimated cutoff" in the output of pseudopotential generation. The current implementation in ld1.x code uses a two-step pseudization: a "hard" norm-conserving pseudopotential is generated, the ultrasoft pseudization is performed on top of it. This yields soft potentials with well-behaving augmentation charges. In order to have a good transferability, at least two channels per angular momentum are required. One channel is typically the bound state; the other is an unbound state, pseudized at some energy in the range of interest for valence electrons. Directories "step_6a" and "step_6b" contain input files, respectively "si_mul.in" and "si_us.in", for the the norm-conserving "hard" pseudopotential and for the ultrasoft pseudopotential.
Note that the two calculations are actually independent: the latter performs the former as well. It is possible to produce norm-conserving, and not only ultrasoft, pseudopotentials, with more than one channel per angular momentum.Look at the difference between the two: apart from the value of "pseudotype", in ultrasoft pseudization you supply a set of matching radii (2.5 in this case), larger than the matching radi for the norm-conserving pseudization (2.0 and 2.1 for s and p respectively). Moreover you need to supply the additional channels for s and p. Try first with -0.1Ry as pseudization energy for both s and p : it will not work, the energy is too far away from the bound state and the algorithm cannot produce a nodeless solution. Try -0.3Ry, without the "rho0=0.01" line : it will not work. By setting a nonzero value of the charge at r=0 (this is what the "rho0" option does) you force a different form of the pseudowave, and this time it works.