A fast and precise DFT wavelet code

Fitting pseudopotentials

A fast and precise DFT wavelet code
Jump to: navigation, search

Note that this tutorial is under construction and subject to changes. All instructions and sample input/output files will be contained in the pseudo/Tutorial/ subdirectory of BigDFT.

This lesson aims to give a short introduction to the use of pseudo, a tool that allows to generate HGH pseudopotentials for use with BigDFT. The step by step instructions address the basic fitting procedure and the inclusion of nonlinear core corrections (NLCC).

Contents

Generating the atomic reference data

  • Before a pseudopotential can be optimized, an atomic all electron DFT calculation is needed. For this we use the tool atom to generate some reference data for the XC functional and atomic configurations of interest.
  • First take a closer look at the following input file atom.dat, which gives two configurations of the atom, one in closed shell, one is spin polarized.
 O
 PADE
spin polarized
 20.0      10.0 60.0             rmax,aa,bb
       1.38d0   6.0d0            rcov,rprb  -> note rprb is rather large, weak confinement
    1    9                       number of core and valence orbitals
    2    0     1.00      1.00   \
    3    0     0.00      0.00    \
    4    0     0.00      0.00    |
    2    1     2.00      2.00    |
    3    1     0.00      0.00     > equal up and down occupations -> closed shell in LSDA
    4    1     0.00      0.00    |
    3    2     0.00      0.00    |
    4    2     0.00      0.00    /
    4    3     0.00      0.00   /
NEXT CONFIGURATION
    2    0     2.00      0.00   \
    3    0     0.00      0.00    \
    4    0     0.00      0.00    |
    2    1     4.00      0.00    |
    3    1     0.00      0.00     > Second channel is zero -> Hund's rule will be used
    4    1     0.00      0.00    |
    3    2     0.00      0.00    |
    4    2     0.00      0.00    /
    4    3     0.00      0.00   /
  • Run the program atom. It will write out some files that include all needed input files for pseudo.
  • Existing Input files for pseudo, such as input.weights and psppar are not overwritten. atom just appends some lines.
  • atom also generates plenty of files for plotting. We will use them later with pseudo and gnuplot.
  • To check whether you got the correct reference data, compare atom.00.ae with atom.00.ae in the ref subdirectory. The file contains eigenvalues and charge integrals up to the covalent radius that will be needed for fitting.
    #    n    l ...          eval                charge
         2    0 ... -0.868941624336E+00  0.725402152717E+00
         2    0 ... -0.868941624336E+00  0.725402152717E+00
         3    0 ...  0.151372211048E-01  0.175632463466E-01
         3    0 ...  0.151372211048E-01  0.175632463466E-01
         4    0 ...  0.883882720213E-01  0.895556274062E-02
         4    0 ...  0.883882720213E-01  0.895556274062E-02
         2    1 ... -0.335890455396E+00  0.652380898401E+00
         2    1 ... -0.335890455396E+00  0.652380898401E+00
         3    1 ...  0.559529537874E-01  0.125193436094E-01
         3    1 ...  0.559529537874E-01  0.125193436094E-01
         4    1 ...  0.112328282101E+00  0.127731262745E-01
         4    1 ...  0.112328282101E+00  0.127731262745E-01
         3    2 ...  0.911350106211E-01  0.649588626749E-04
         3    2 ...  0.911350106211E-01  0.649588626749E-04
         4    2 ...  0.138725530583E+00  0.183082913971E-03
         4    2 ...  0.138725530583E+00  0.183082913971E-03
         4    3 ...  0.122461690116E+00  0.220707980970E-06
         4    3 ...  0.122461690116E+00  0.220707980970E-06

Testing the input guess pseudopotential

  • Take a look at input.pseudo generated by the program atom. Verify that the kewyord -fit is not present, as we want to run pseudo first without fitting.
  • Copy psppar.GTH to psppar and run pseudo in serial to test the accuracy of this pseudopotential for a closed shell atom. You can find this pseudopotential and many others on the download pages of abinit].
this is the Goedecker-Teter-Hutter pseudopotential for Oxygen
8   6   960531                      zatom,zion,pspdat
2   1   0   0  2001  0.             pspcod,pspxc,lmax,lloc,mmax,r2well
0.2477535 -16.4822284 2.3701353 0 0 rloc, c1, c2, c3, c4
.2222028 18.1996387 0               rs, h1s, h2s
0 0                                 rp, h1p
  • Compare your output with pseudo.out.ref in the ref subdirectory - it should be identical.
  Pseudo atom energies
  ____________________

 kinetic energy            =   11.8476868721
 potential energy          =  -16.2328629373
 non local energy          =    1.3038686713
 hartree energy with XC    =   11.3473668807
 exchange + corr energy    =   -3.1536086331
 vxc    correction         =   -4.1347595852
 -------------------------------------------
 el-el  interaction energy =   13.6323865041
 external energy           =  -39.3628763547
 sum of eigenvalues        =   -3.0813073939
 -------------------------------------------
 total energy              =  -15.7325429459

Side note: The accuracies are relatively poor here because we have a weaker confinement (larger rprb) than the one used for the GTH oxygen pseudopotential.

  • One should always take a look at the orbital plots. You can do this with gnuplot using the script orbitals.gplt.
gnuplot orbitals.gplt

Running a first short optimization

  • Add the flags -fit and -maxiter 100 to the file input.pseudo, which results in a short optimization of the psppar.
  • -maxiter is only used to keep the output short here, usually no such option is needed.
  • Run pseudo in serial and observe how the penalty function is decreasing.
  • The output should look very similar to pseudo.shortfit.out (with deviations from random seeds ...)
 nl    s      occ              ae        pseudo       diff     diff*weight
 2s  -0.5    1.0000
         eigenvalue      -0.8689E+00 -0.8689E+00  0.1018E-07  0.1018E-03
         charge           0.7254E+00  0.7254E+00  0.1262E-05  0.1262E-02
         residue                      0.2271E-07              0.0000E+00
         rnode                        0.0000E+00              0.0000E+00
  • The actual choice of the penalty function in input.weights can have a big effect on the result of your optimization. For now we will use the given example.
0d0 0d5 0d0 0d0 0d0 1d3  psi(0), dEkin_wvlt, radii, hij, locality  E_exct
   n   l  so    eigval  chrg    dchrg  ddchrg  res    rnode dnode ddnode
   2   0 -0.50   1.0e4   1.0e3   0.0e0  0.0e0   0.0e0  1.0e0 0.0e0 0.0e0
   2   0  0.50   1.0e4   1.0e3   0.0e0  0.0e0   0.0e0  1.0e0 0.0e0 0.0e0
   3   0 -0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   3   0  0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   0 -0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   0  0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   2   1 -0.50   1.0e4   1.0e3   0.0e0  0.0e0   0.0e0  1.0e0 0.0e0 0.0e0
   2   1  0.50   1.0e4   1.0e3   0.0e0  0.0e0   0.0e0  1.0e0 0.0e0 0.0e0
   3   1 -0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   3   1  0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   1 -0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   1  0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   3   2 -0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   3   2  0.50   1.0e0   1.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   2 -0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   2  0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   3 -0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0
   4   3  0.50   0.0e0   0.0e0   0.0e0  0.0e0   0.0e0  0.0e0 0.0e0 0.0e0

Invoking pseudo in parallel

  • Remove the -fit flag from input.pseudo to do another test without fitting.
  • Run pseudo in parallel using two MPI processes, for example
mpirun -np 2 $CompileBigDFT/pseudo/src/pseudo | tee test.out

One process will treat the configuration atom.00.ae, the other one atom.01.ae.

  • The standard output should give the same eigenvalues, charges, etc. as it did before.
  • The file proc.01.out gives the corresponding output for the polarized configuration atom.01.ae, like ref/4.proc.01.out.
 nl    s      occ              ae        pseudo       diff     diff*weight
 2s  -0.5    1.0000
         eigenvalue      -0.9135E+00 -0.9170E+00 -0.3452E-02  0.3452E+02
         charge           0.7333E+00  0.7348E+00  0.1529E-02  0.1529E+01
         residue                      0.2794E-07              0.0000E+00
         rnode                        0.0000E+00              0.0000E+00

 2s   0.5    1.0000
         eigenvalue      -0.7987E+00 -0.7928E+00  0.5923E-02  0.5923E+02
         charge           0.7112E+00  0.7083E+00 -0.2859E-02  0.2859E+01
         residue                      0.1658E-07              0.0000E+00
         rnode                       -0.6790E-04              0.6790E-04
...
 2p  -0.5    3.0000
         eigenvalue      -0.3791E+00 -0.3829E+00 -0.3880E-02  0.3880E+02
         charge           0.6706E+00  0.6739E+00  0.3222E-02  0.3222E+01
         residue                      0.7274E-06              0.0000E+00
         rnode                        0.0000E+00              0.0000E+00

 2p   0.5    1.0000
         eigenvalue      -0.2694E+00 -0.2626E+00  0.6748E-02  0.6748E+02
         charge           0.6181E+00  0.6108E+00 -0.7278E-02  0.7278E+01
         residue                      0.3001E-06              0.0000E+00
         rnode                        0.0000E+00              0.0000E+00
  • Note that pseudo also reports energy differences of the configurations, which in this case is the atomic spin polarization energy. We will need NLCC later to describe spin effects like this quantity.
 Excitation energies
 ___________________

 Configuration,    dE=Etot-Etot1,    (dE-dE_AE)*weight
             0  0.000000000000E+00  0.000000000000E+00
             1 -0.605168248675E-01  0.543678408899E+01

You can fit the energy differences with the weight E_exct in input.weights.

0d0 0d5 0d0 0d0 0d0 1d3  psi(0), dEkin_wvlt, radii, hij, locality  E_exct

Adding NLCC to a pseudopotential

  • Compare your pseudopotential with psppar.addNLCC.
spin-polarized     1.38        6.00        method, rcov and rprb
    8   6  20130406                        zatom, zion, date (yymmdd)
   11       -20 0 0 2002 0                 pspcod, IXC, lmax, lloc, mmax, r2well
        0.24775350  2    -16.50052797      2.37488057      rloc nloc c1 .. cnloc
    1                                      nsep
        0.22217380  1     18.20195442      s-projector
        0.14               0.45

Note that pspcod on line 3 takes a value of 11, which means NLCC are used.

  • The last line gives a first input guess for the width rcore and the charge fraction qcore for the NLCC, respectively. Update your psppar accordingly.
  • Run pseudo without fitting to see the resulting orbital properties. The pseudopotential is not optimized for the NLCC, so the result is quite bad.
  • Note that a file nlcc.gnuplt is written that allows you to plot and fit the core charge for the NLCC with gnuplot.
gnuplot
load 'nlcc.gplt'
  • Run pseudo with -fit and -maxiter 200 in input.pseudo using two MPI processes. Compare your result with ref/pseudo.NLCCfit.out.ref.
 Excitation energies
 ___________________

 Configuration,    dE=Etot-Etot1,    (dE-dE_AE)*weight
             1 -0.567411526129E-01  0.166111183439E+01

Customizing the free fitting parameters

  • When pseudo reads input.fitpar, it reports which parameters are fitted. These names can be copied to input.fitpar to free certain parameters.
  • The keyword auto in input.fitpar frees all fitting parameters except for rloc, rcore and qcore.
  • Add the keywords rcore and qcore on any line in input.fitpar such that the NLCC can be fitted as well.
 Reading fitting parameters from input.fitpar
 ____________________________________________

 auto: Free all params except for rloc, rcore and qcore
 qcore corresponds to gcore1
  2 lines have been read. Free params are:
 gpot1
 gpot2
 rcore
 gcore1
 rs
 hs11
  • Do another fit using two MPI processes without the -maxiter keyword. The penalty for the excitation energies should improve quite a bit.
 Excitation energies
 ___________________

 Configuration,    dE=Etot-Etot1,    (dE-dE_AE)*weight
             1 -0.550568319213E-01  0.232088572076E-01
  • Plot the core charges again using nlcc.gnuplt and compare your results with 6.fitNLCC.out.
Personal tools