Home > Scripts > Models > Factory > sqw_phonons.m



model=sqw_phonons(configuration, calculator, ..., options)


function signal=sqw_phonons(configuration, varargin)


 model=sqw_phonons(configuration, calculator, ..., options)

   iFunc/sqw_phonons: computes phonon dispersions using the ASE.
   A model which computes phonon dispersions from the Hellmann-Feynman forces acting 
     between atoms. The input argument is any configuration file describing the
     material, e.g. CIF, PDB, POSCAR, ... supported by ASE. Alternatively, a
     Crystallography Open Database ID or chemical formulae can be entered.

   This models can compute the coherent inelastic phonons dispersions for
     any crystalline (powder or single crystal) material, in the harmonic and
     adiabatic approxiimation.
   The phonon spectra is computed using one of the calculator supported by the
   Atomic Simulation Environment (ASE) <https://wiki.fysik.dtu.dk/ase>.

   Supported calculators are (use name as input argument, or set 'calculator=...'):
     ABINIT    Plane-wave pseudopotential code
     ELK       Full Potential LAPW code
     EMT       Effective Medium Theory calculator (Al,Cu,Ag,Au,Ni,Pd,Pt,H,C,N,O)
     GPAW      Real-space/plane-wave/LCAO PAW code
     QuantumEspresso Plane-wave pseudopotential code
     VASP      Plane-wave PAW code (when installed and licensed)

   We recommend QuantumEspresso, GPAW, ABINIT and Elk.
   Refer to https://molmod.ugent.be/deltacodesdft for codes comparison.

   The simplest usage is to call: sqw_phonons('gui') which displays an entry dialog box
   and proceeds with a fully automatic computation, and plots final results.
   The calculators can be specified by just giving their name as a parameter, 
   or using e.g. options.calculator='GPAW'. Except for EMT, other calculators must
   be installed separately, and optionally specific pseudo-potentials. 
     See https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html

   Benchmarks indicate that, for phonon dispersions:
   * QuantumEspresso, with excellent parallelization and accuracy.
   * ABINIT, VASP are also fast codes.
   * The all-electrons Elk code is about 10 times slower than QuantumEspresso.
   * Using k-points=3 is both fast and ensures a reasonable accuracy in phonons
   * Phonon dispersions are not too sensitive on the energy cut-off. 340 eV is good.
   * Elk and ABINIT can not handle large supercells without recompiling.
   * NWChem and Elk are not sensitive to the energy cut-off.

 The model must first be created (which triggers e.g. a DFT computation), and 
 once created, it can be evaluated in the whole HKLE reciprocal space.

 MODEL CREATION ===============================================================

 sqw_phonons(configuration, calculator, options, ...) 
   is the usual syntax, such as sqw_phonons('POSCAR','qe','metal')
   shows a dialogue to select the phonon calculator configuration.
 sqw_phonons('defaults') or sqw_phonons
   uses a simple fcc Aluminium as example.
 sqw_phonons(dir, calculator, options, ...)
   re-use an existing directory containing any of POSCAR, ASE atom object, PhonoPy
   files, ...

 The arguments for the model creation should be:

 input (model creation):
 configuration: file name or directory to an existing material configuration
   Any A.S.E supported format can be used (POSCAR, CIF, SHELX, PDB, ...). 
     See <https://wiki.fysik.dtu.dk/ase/ase/io.html#module-ase.io>
   Alternatively, the 'bulk','molecule', 'crystal', and 'nanotube' ASE constructors can be
   used, using the Python syntax, e.g. 
       'bulk("Si", "diamond", a=5.4)'
       'bulk("Cu", "fcc", a=3.6, cubic=True)'
       'nanotube(6, 0, length=4)'
       'crystal(["Na", "Cl"], [(0, 0, 0), (0.5, 0.5, 0.5)], spacegroup=225,
          cellpar=[5.64, 5.64, 5.64, 90, 90, 90])'
     See <https://wiki.fysik.dtu.dk/ase/ase/structure.html>
   Alternatively, the structure can be entered as a COD ID such as
       'cod: 1532356' will fetch ID on Crystallography Open Database.
   Alternatively, the structure can be entered as a chemical formulae in Hill
   notation (C, H and then alpha order) on Crystallography Open Database. When
   multiple matches, a dialogue will pop-up to select structure.
       'cod: O7 Se2 Tb2'

 'metal','insulator','semiconductor': indicates the type of occupation for 
    electronic states, which sets smearing.

 'html' or 'report':   generate a full report of the simulation and export data.
 'plot' or 'autoplot': plot the results after building the model.

 options: an optional structure with optional settings, as follows:

 General options
   options.target =path                   Where to store all files and FORCES
     a temporary directory is created when not specified.
   options.supercell=scalar or [nx ny nz] supercell size. Default is 0 (auto mode).
   options.calculator=string              EMT,GPAW,Elk,ABINIT,Quantum
     We recommend ABINIT,QE and Elk. Default set from installed software.
     The calculator can also be given as a Python string, which then 
     also requires to import the proper Python module. Use e.g.:
       options.declaration = 'from ase.calculators.lammpsrun import LAMMPS';
       options.calculator  = 'LAMMPS()';
       options.command     = '/usr/bin/lammps';
   options.potentials=string              Basis datasets or pseudopotentials.
     Elk:    path to potentials, e.g. /usr/share/elk-lapw/species (ELK_SPECIES_PATH)
     ABINIT: path to potentials, e.g. /usr/share/abinit/psp/      (ABINIT_PP_PATH)
       or any of 'NC' 'fhi' 'hgh' 'hgh.sc' 'hgh.k' 'tm' 'paw'
     QuantumEspresso: path to potentials, e.g. /usr/share/espresso/pseudo
   options.command='exe'                  Path to calculator executable.
   options.mpi=scalar                     use multi-cpu, requires e.g. OpenMPI to be installed.
     when available and not specified, MPI is used with all CPU's.
   options.machinefile=filename           file containing the list of MPI machines to use
   options.autoplot=0|1                   when set, automatically evaluates the Model
     and plot results.
   options.htmlreport=0|1                 when set, automatically generates a full
     report on the computation results.
   options.optimizer                      when set, performs a geometry optimization
     before computing the forces (move atoms in the cell). Can be set to 
       'BFGS' or 'LBFGS' (low memory BFGS)
       'MDMin' or 'FIRE'
   options.accuracy                       'fast', 'very fast' (default) or 'accurate'.
     The 'fast' choice uses the symmetry operators to lower the number of atom 
     displacements. The force gradient uses central difference.
     The 'very fast' option halves the number of displacements. The equilibrium
     forces are used to improve the accuracy of the force gradient.
     When using the 'accurate','fast' or 'very fast' accuracy, the calculation
     is limited to 256 iteration steps. For more, use 'single_shot' mode below.
     The 'single_shot' option performs all in a single call. No ETA is displayed, 
     but computation is the fastest, and allows more than 256 iterations.
     The 'accurate' choice is longer to execute (e.g. 3-6 times slower), but computes all forces.
   options.use_phonopy=0|1                requests to use PhonoPy when installed.
     This choice also sets accuracy='very fast'
   options.disp=value                     the atom displacement in Angs. Default is 0.01.

 DFT specific options
   options.kpoints=scalar or [nx ny nz]   Monkhorst-Pack grid. Default is 0 (auto mode).
   options.xc=string                      Type of Exchange-Correlation functional to use
     'LDA','PBE','revPBE','RPBE','PBE0','B3LYP'   for GPAW
     'LDA','PBE','revPBE','RPBE','GGA'            for ABINIT
     'LDA','PBE','REVPBE','PBESOL','WC06','AM05'  for ELK
     'LDA','PBE','PW91'                           for VASP
     Default is 'PBE'.
   options.raw='name=value, ...'          Additional arguments passed to the ASE
                                          using the Python syntax.
   options.nbands=int                     Number of valence bands for which 
     wavefunctions are being computed. (= # of electrons /2); for a metal, 20% more
     (minimum 4 more).
   options.ecut=scalar                    Kinetic energy cutoff for wavefunctions. 
     Large value improves convergence. Estimate is 200*n_typ_atoms in [eV].
     Usually one runs at calculations at various ecut to investigate convergence
     ABINIT paw/pawxml potentials require larger ecut (1000 eV) and nsteps (100).
   options.diagonalization='dav' or 'cg' or 'rmm-diis'
     The Davidson method is faster than the conjugate-gradient but uses more
     memory. The RMM-DIIS method allows parallelization.
   options.occupations='metal'            for metals ('smearing') help converge
                       'insulator'        for insulators
                       'semiconductor'    sets 0 eV FermiDirac smearing
                       'auto'             for Elk (automatic smearing)
                       or 0 for semi-conductors
                       or a value in eV for the distribution width e.g. 0.3
   options.nsteps=scalar                   Max number of iterations for SCF.
     Typical: 30. Large value improves convergence.
   options.toldfe=scalar                  Convergence threshold on the energy for 
     selfconsistency, e.g. 1e-7 [eV]. Default is 1e-5 [eV] for fast estimate.

 Options specific per calculator
   options.mode='pw','fd', or 'lcao'      GPAW computation mode as Plane-Wave,
     Finite Difference, or LCAO (linear combination of atomic orbitals). Default is 'pw'.
   options.iscf='NC','PAW'                Type of SCF cycles (ABINIT) 
   options.pps = 'fhi' 'hgh' 'hgh.sc' 'hgh.k' 'tm' 'paw' 'pawxml'      (ABINIT)
                 a structure such as struct('Ca','_pv') ...              (VASP)
                 'fhi' 'tm'                                            (SIESTA)
                 'hscv_pbe' 'hscv_lda' 'hgh_lda' 'sg15'               (Octopus)
   options.mixing_beta=scalar             mixing factor for self-consistency
     default=0.7. use 0.3 to improve convergence (QuantumEspresso)
   options.mixing_ndim=scalar             number of iterations used in mixing
     default=8. If you are tight with memory, you may reduce it to 4. A larger
     value will improve the SCF convergence, and use more memory (QuantumEspresso)

 Options to specify executables
   options.calculator                     Python/ASE calculator (see above)
   options.command='exe'                  Path to calculator executable

   and executables can be specified for 'python','mpirun','gpaw','elk','abinit',
     'quantumespresso','vasp','octopus','cp2k','siesta', for instance:

 The options can also be entered as strings with 'field=value; ...'.

 output (model creation):
   model: iFunc 4D model S(qh,qk,ql,w) with parameters p.

 The syntax sqw_phonons(model,'html') allows to re-create the HTML report about
 the 4D phonon model.

 You may look at the following resources to get material structure files:

   This model is suitable to compute phonon dispersions for e.g solid-
   state materials.
   The Atomic Simulation Environment must be installed.
   The temporary directories (UserData.dir) are not removed.
   Momentum is expressed in reduced lattice unit, and energy is in meV.
        1 meV = 241.8 GHz = 11.604 K = 0.0965 kJ/mol = 8.0657 cm-1 

 References: https://en.wikipedia.org/wiki/Phonon

 Atomic Simulation Environment
   S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng., Vol. 4, 56-66, 2002
   <https://wiki.fysik.dtu.dk/ase>. Exists as Debian package 'python-ase'.
   Repository: http://download.opensuse.org/repositories/home:/dtufys/
 GPAW J. J. Mortensen et al, Phys Rev B, Vol. 71, 035109 (2005). 
   <http://wiki.fysik.dtu.dk/gpaw>. Exists as Debian package 'gpaw' and 'gpaw-setups' (gpaw-data).
 Elk <http://elk.sourceforge.net>. Exists as Debian package 'elk-lapw'. 
   The Elk executable should be compiled with e.g. 
     elk/src/modmain.f90:289 maxsymcrys=1024 or larger
 ABINIT  X. Gonze et al, Computer Physics Communications 180, 2582-2615 (2009).
   <http://www.abinit.org/>. Exists as Debian package 'abinit' and 'abinit-doc' (abinit-data).
   Install potentials from http://wiki.fysik.dtu.dk/abinit-files/abinit-pseudopotentials-2.tar.gz
   into e.g. /usr/share/abinit
 PHON D. Alfe, Computer Physics Communications 180,2622-2633 (2009). 
 Quantum Espresso P. Giannozzi, et al, J.Phys.:Condens.Matter, 21, 395502 (2009).
 VASP G. Kresse and J. Hafner. Phys. Rev. B, 47:558, 1993.
   <https://www.vasp.at/> Requires a license.
 PhonoPy:   A. Togo and I. Tanaka, Scr. Mater., 108, 1-5 (2015)

 MODEL EVALUATION (once created) ==============================================

 Once the model has been created, its use requires that axes are given on
   regular qx,qy,qz grids (in rlu along reciprocal axes). The model evaluations
   does not require to recompute the forces, and is very fast. 
 When all axes are vectors of same orientation, the HKL locations is assumed to be a q-path.
 When axes are not all vectors, not same length, nor orientation, a 3D HKL cube is used.
 To generate a powder 2D S(q,w) you may use: powder(model)
 To generate the dispersion curves along principal crystallographic directions,
   use: band_structure(model)

 The phonon density of states (DOS) is obtained with: dos(model)
 by evaluation onto a grid HKL=[-.5 : 0.5]. It is stored in model.UserData.DOS as an 
 iData object. Partial density of states (per mode) are stored in UserData.DOS_partials

 The neutron scattering cross sections should be automatically determined from 
 the chemical formula in order to properly evaluate the neutron scattering 
 In case this would fail, or other cross sections are needed, it is required
 to set the coherent neutron scattering length (b_coh in [fm]) as a vector in the 
   model.UserData.properties.b_coh = [vector / number of atoms in the cell]
 or alternatively the scattering cross section in [barn]
   model.UserData.properties.sigma_coh = [vector / number of atoms in the cell]
 where negative values set b_coh < 0 (e.g. Hydrogen).
 Setting b_coh=0 unactivates the intensity estimate.
 This assignement should be done after creating the model, before performing
 further HKL evaluations. Default is to use guessed b_coh from the formula.

 Once the model has been created:
 input:  p: sqw_phonons model parameters (double)
             p(2)=Gamma   dispersion DHO half-width in energy [meV]
             p(3)=Background (constant)
             p(4)=Temperature of the material [K]. When 0, the intensity is not computed.
             p(5)=Energy_scaling. All frequencies are multiplied by this value.
          or p='guess'
         qh: axis along QH in rlu (row,double)
         qk: axis along QK in rlu (column,double)
         ql: axis along QL in rlu (page,double)
         w:  axis along energy in meV (double)
    signal: when values are given, a guess of the parameters is performed (double)
 output: signal: model value

 Example (model creation and evaluation):
   s=sqw_phonons('bulk("Cu", "fcc", a=3.6, cubic=True)','EMT','metal');
   plot3(s)            QH=0 plot of the Brillouin Zone

   publish(s)          generates a full report for the given model.
   band_structure(s)   plots the dispersion curves and density of states.

   s=sqw_phonons('bulk("Si", "diamond", a=5.4, cubic=True)','semiconductor');

   s=sqw_phonons([ ifitpath 'Data/POSCAR_Al'],'metal','EMT');

 Version: Nov. 26, 2018
 See also iData, iFunc/fits, iFunc/plot, gauss, sqw_cubic_monoatomic, sqw_sine3d, sqw_vaks
   iFunc_Sqw4D/powder, <a href="matlab:doc(iFunc,'Models')">iFunc:Models</a>
 (c) E.Farhi, ILL. License: EUPL.


This function calls: This function is called by:
Generated on Mon 26-Nov-2018 15:08:42 by m2html © 2005. iFit (c) E.Farhi/ILL EUPL 1.1