A fast and precise DFT wavelet code

First runs with BigDFT

A fast and precise DFT wavelet code
Revision as of 14:34, 11 June 2018 by Deutsch (Talk | contribs)

Jump to: navigation, search

This lesson has been created for the current stable version. Earlier versions are fully capable of running this tutorial but input files may have to be changed according to possible earlier formats.

Basics of BigDFT: first runs and managing different calculations, N2 molecule as example

BigDFT code is organized by optional input files, each one with compulsory variables. Input files all have a specific extension (.xyz, .geopt, .kpt ...), each one being associated to a particular type of calculation (The latest version of BigDFT now also accepts .yaml files in which you can directly put all the informations needed regarding the run you want to perform). With this lesson you will have to deal with the different outputs of BigDFT code, such as to learn how to manipulate basic DFT objects.

Default run: no input files

As mentioned above, all the BigDFT files are optional, except the file of the atomic positions. This means that a DFT calculation can be done also by giving the atomic positions file only. Consider for example the N2 molecule, given by the posinp.xyz file:

2 angstroem
N 0. 0. 0.
N 0. 0. 1.11499

Run the code in a directory which has only this file.

user@garulfo:~/N2/$ ls
bigdft  posinp.xyz
user@garulfo:~/N2/$ ./bigdft | tee N2.out

The screen output should then behave like that:

  #------------------------------------------------------------------------ Input parameters
   debug                               : No #     debug option
   fftcache                            : 8192 #   cache size for the FFT
   accel                               : NO #     acceleration
   ocl_platform                        : ~ #      Chosen OCL platform
   ocl_devices                         : ~ #      Chosen OCL devices
   blas                                : No #     CUBLAS acceleration
   projrad                             : 15. #    Radius of the projector as a function of the maxrad
   exctxpar                            : OP2P #   Exact exchange parallelisation scheme
   ig_diag                             : Yes #    Input guess (T=Direct, F=Iterative) diag. of Ham.
   ig_norbp                            : 5 #      Input guess Orbitals per process for iterative diag.
   ig_blocks: [300, 800] #                        Input guess Block sizes for orthonormalisation
   ig_tol                              : 1e-4 #   Input guess Tolerance criterion
   methortho                           : 0 #      Orthogonalisation
   rho_commun                          : DEF #    Density communication scheme (DBL, RSC, MIX)
   psolver_groupsize                   : 0 #      Size of Poisson Solver taskgroups (0=nproc)
   psolver_accel                       : 0 #      Acceleration of the Poisson Solver (0=none, 1=CUDA)
   unblock_comms                       : OFF #    Overlap Communications of fields (OFF,DEN,POT)
   linear                              : OFF #    Linear Input Guess approach
   tolsym                              : 1e-8 #   Tolerance for symmetry detection
   signaling                           : No #     Expose calculation results on Network
   signaltimeout                       : 0 #      Time out on startup for signal connection (in seconds)
   domain                              : ~ #      Domain to add to the hostname to find the IP
   inguess_geopt                       : 0 #      input guess to be used during the optimization
   store_index                         : Yes #    store indices or recalculate them for linear scaling
   verbosity                           : 2 #      verbosity of the output
   outdir                              : . #      Writing directory
   psp_onfly                           : Yes #    Calculate pseudopotential projectors on the fly
   pdsyev_blocksize                    : -8 #     SCALAPACK linear scaling blocksize
   pdgemm_blocksize                    : -8 #     SCALAPACK linear scaling blocksize
   maxproc_pdsyev                      : 4 #      SCALAPACK linear scaling max num procs
   maxproc_pdgemm                      : 4 #      SCALAPACK linear scaling max num procs
   ef_interpol_det                     : 1e-20 #  FOE max determinant of cubic interpolation matrix
   ef_interpol_chargediff              : 10. #    FOE max charge difference for interpolation
   mixing_after_inputguess             : Yes #    mixing step after linear input guess (T/F)
   iterative_orthogonalization         : No #     iterative_orthogonalization for input guess orbitals
   check_sumrho                        : 2 #      enables linear sumrho check
   experimental_mode                   : No #     activate the experimental mode in linear scaling
   write_orbitals                      : No #     linear scaling write KS orbitals for cubic restart (might take lot of disk space!)
   hgrids: [0.45, 0.45, 0.45] #                   grid spacing in the three directions (bohr)
   rmult: [5., 8.] #                              c(f)rmult*radii_cf(:,1(2))=coarse(fine) atom-based radius
   ixc                                 : 1 #      exchange-correlation parameter (LDA=1,PBE=11)
   ncharge                             : 0 #      charge of the system
   elecfield: [0., 0., 0.] #                      electric field (Ex,Ey,Ez)
   nspin                               : 1 #      spin polarization
   mpol                                : 0 #      total magnetic moment
   gnrm_cv                             : 1e-4 #   convergence criterion gradient
   itermax                             : 50 #     max.
   nrepmax                             : 1 #      max.
   ncong                               : 6
   idsx                                : 6 #      wfn. diis history
   dispersion                          : 0 #      dispersion correction potential (values 1,2,3,4,5), 0=none
   inputpsiid                          : 0
   output_wf                           : 0
   output_denspot                      : 0
   rbuf                                : 0. #     length of the tail (AU)
   ncongt                              : 30
   norbv                               : 0 #      Davidson subspace dim.
   nvirt                               : 0
   nplot                               : 0
   disablesym                          : No #     disable the symmetry detection
   method                              : manual # K-point sampling method
   kpt: #                                         Kpt coordinates
   -  [0., 0., 0.]
   wkpt: [1.] #                                   Kpt weights
   bands                               : No #     For doing band structure calculation
   method                              : none #   Geometry optimisation method
   ncount_cluster_x                    : 1 #      Maximum number of force evaluations
   frac_fluct                          : 1.
   forcemax                            : 0.
   randdis                             : 0. #     random displacement amplitude
   betax                               : 4. #     Stepsize for the geometry optimisation
   iscf                                : 0 #      mixing parameters
   itrpmax                             : 1 #      maximum number of diagonalisation iterations
   rpnrm_cv                            : 1e-4 #   stop criterion on the residue of potential or density
   norbsempty                          : 0 #      No. of additional bands
   Tel                                 : 0. #     electronic temperature
   occopt                              : 1 #      smearing method
   alphamix                            : 0. #     Multiplying factors for the mixing
   alphadiis                           : 2. #     Multiplying factors for the electronic DIIS
   sic_approach                        : none #   SIC method
   sic_alpha                           : 0. #     SIC downscaling parameter
   tddft_approach                      : none #   TDDFT method
  #--------------------------------------------------------------------------------------- |
 Data Writing directory                : ./
  #-------------------------------------------------- Input Atomic System (file: posinp.xyz)
 Atomic System Properties:
   Number of atomic types              :  1
   Number of atoms                     :  2
   Types of atoms                      :  [ N ]
   Boundary Conditions                 : Free #Code: F
   Number of Symmetries                :  0
   Space group                         : disabled
  #------------------------------ Geometry optimization Input Parameters (file: input.geopt)
 Geometry Optimization Parameters:
   Maximum steps                       :  1
   Algorithm                           : none
   Random atomic displacement          :  0.0E+00
   Fluctuation in forces               :  1.0E+00
   Maximum in forces                   :  0.0E+00
   Steepest descent step               :  4.0E+00
 Material acceleration                 :  No #iproc=0

In the log file, BigDFT automatically displays all the input parameters used for the calculations (this is what you see above). The parameters which were not explicitly given are set to a default value, except the atomic positions of course, which has to be given. Each input file contains compulsory lines, with the exception of input.perf which controls developer-oriented performance variables. You can see there the possible optional files which BigDFT might read . Since they do not exist, their default values are applied to the code. Basically, they correspond to a single-point LDA calculation, without k-points nor spin-polarisation. As we said earlier, the new version of BigDFT now supports .yaml files were all the parameters are given (instead of using an optional file every time you want to calculate something specific, you can now modify the corresponding line in the .yaml file), the input display shown above can serve as a template for this .yaml file. -note that one can combine a .yaml file with other task specifying files (.xyz or .geopt for example), BigDFT will still understand-

Using a naming scheme for IO files

All input parameters can be found in files with a naming prefix. By default, this prefix is input (or posinp for atomic input positions). For instance, parameters for geometry optimization will be set up by a file named input.geopt. One can choose the naming prefix by providing an argument to bigdft command line.

Upon launching, BigDFT will look by default for those input (and posinp) files to see which parameters it should use, but one can also decide to rename all those files with a particular name for more practicity. One then has to give this name when launching, so that BigDFT knows what to look for.

Imagine for example that you are interested in visualizing the wavefunctions output of the calculation. To do that, you should enter the suitable parameters in the .dft file (or .yaml file). Let us do it : Create a new calculation set by using the "LDA" prefix and rename all relevant files with LDA:

user@garulfo:~/N2/$ cp posinp.xyz LDA.xyz

Create LDA.dft by copying/pasting in a file named LDA.dft the few lines that follow :

0.550 0.550 0.550   hx,hy,hz: grid spacing in the three directions
3.5  9.0            crmult, frmult: c(f)rmult*radii_cf(*,1(2)) gives the coarse (fine)radius around each atom
1                   ixc: exchange-correlation parameter (LDA=1,PBE=11)
0 0.0  0.0  0.0     ncharge: charge of the system, Electric field
1  0                nspin=1 non-spin polarization, mpol=total magnetic moment
1.E-04              gnrm_cv: convergence criterion gradient
50 10               itermax,nrepmax: maximum number of wavefunction optimizations and of re-diagonalised runs
6  6                ncong, idsx: # CG iterations for the preconditioning equation, length of the diis history
0                   dispersion correction functional (values 1,2,3), 0=no correction
0 0   0             InputPsiId, output_wf, output_denspot
0.0  30             rbuf, ncongt: length of the tail (AU),# tail CG iterations
0 0 0               davidson treatment, no. of virtual orbitals, no of plotted orbitals
T                   disable the symmetry detection

Now you have your default file! It is time to modify LDA.dft such as to output the wavefunctions at the end of calculation, by putting the output_wf variable to 1

(Which repsesents the formatted output. Also 2 (binary) and 3 (ETSF format) are available)

 0 1 0          InputPsiId, output_wf, output_grid 

Now you can run this input file, by putting "LDA" as a command line argument of the code:

user@garulfo:~/N2/$ ./bigdft LDA | tee LDA.out

You can now see that the LDA.dft file is read.

When using a naming scheme, the output files are placed in a directory called data-{naming scheme}. In our LDA example, the wavefunctions of the system can thus be found in the data-LDA directory:

user@garulfo:$ ls data-LDA/
wavefunction-k001-NR.b0001  wavefunction-k001-NR.b0002  wavefunction-k001-NR.b0003
wavefunction-k001-NR.b0004  wavefunction-k001-NR.b0005

Here 001 means the first K-point (meaningless in this case), N stands for non spin-polarized, R for real part and the remaining number is the orbital ID. Post-processing of these files is done in the fourth tutorial.

In the same spirit, another calculation can be done with different parameters. Imagine we want to perform a Hartree-Fock calculation. In BigDFT, this can be done by putting the ixc input variable to 100. So, copy the LDA.dft file to HF.dft and modify it accordingly (don't forget to rename also the coordinate file). This time, an error will occur:

user@garulfo:$ ./bigdft HF | tee HF.out
 ERROR: The pseudopotential parameter file "psppar.N" is lacking,
        and no registered pseudo found for "N", exiting...

This is because the pseudopotential is assigned by default in the code only for LDA and PBE XC approximations. You can find here the pseudopotential which is taken by default in the LDA run. Put it in a psppar.N file, and run the calculation (redirect the output in the file HF.out for example). When possible, care should be taken in choosing a pseudopotential which has been generated with the same XC approximation used. Unfortunately, at present HGH data are only available for semilocal functionals. For example, the same exercise as follows could have been done with Hybrid XC functionals, like for example PBE0 (ixc=-406). See the XC codes for a list of the supported XC functionals. Data of the calculation can be analysed. Consider the eigenvectors of the LDA and HF runs.

Exercise: Compare the values of the HOMO and HOMO-1 eigenvalues for the LDA and the HF run. Change the values of the hgrid and crmult to find the converged values. Note that, both in the LDA and in the HF calculation, a norm-conserving PSP is used. The results can be compared to all-electron calculations, done with different basis sets, from references (units are eV) [1] and [2]:

        LDA(1)      HF(1)      HF(2)       (Exp.)
3σg    10.36        17.25      17.31      (15.60) 
1πu    11.84        16.71      17.02      (16.98)
2σu    13.41        21.25      21.08      (18.78)

The results depends, of course, on the precision chosen for the calculation, and of the presence of the pseudopotential. As it is well-known, the pseudopotential appoximation is however much less severe than the approximation induced by typical XC functionals. We might see that, even in the HF case, the presence of a LDA-based pseudopotential (of rather good quality) does not alter so much the results. Here you can find the values from BigDFT calculation using a very good precision (hgrid=0.3, crmult=7.0). Note that 1 ha=27.21138386 eV.

        LDA           HF     
3σg    10.40        17.32     
1πu    11.75        16.62    
2σu    13.52        21.30

How much do these values differ from the calculation with default parameters? Do they converge to a given value? What is the correlation for the N2 molecule in (PSP) LDA?

NOTE: when you create your own .xyz files, be sure to write the names of the atoms in capital letters .

  1. S. Hamel et al. J. Electron Spectrospcopy and Related Phenomena 123 (2002) 345-363
  2. P. Politzer, F. Abu-Awwad, Theor. Chem. Acc. (1998), 99:83-87
Personal tools