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.
   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:
     ABINIT    Plane-wave pseudopotential code
     Dacapo    Plane-wave ultra-soft 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
     NWChem    Gaussian based electronic structure code
     QuantumEspresso_ASE Plane-wave pseudopotential code (with ASE/QE-util)
     QuantumEspresso     Plane-wave pseudopotential code (with PHON)
     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/PHON is the fastest, with excellent parallelization and accuracy.
   * QuantumEspresso/ASE is excellent.
   * 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 ===============================================================

 The arguments for the model creation should be:

 input (model creation):
 configuration: file name 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>

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

 'dos': will also compute the vibrational density of states.
 '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

 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,NWChem,Dacapo,ABINIT,Quantum
     We recommend ABINIT,QE and Elk. Default set from installed software.
   options.dos=1                          Option to compute the vibrational
     density of states (vDOS) in UserData.DOS
   options.potentials=string              Basis datasets or pseudopotentials.
     NWChem: see http://www.nwchem-sw.org/index.php/Release64:AvailableBasisSets
       e.g. 'ANO-RCC'
     Dacapo: path to potentials, e.g. /usr/share/dacapo-psp       (DACAPOPATH)
     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.
   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. Also requests vDOS computation (options.dos=1).
   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                       can be 'fast' (default) or 'accurate'.
     the 'fast' choice uses the symmetry operators to lower the number of atom 
     displacements. This lowers the accuracy of the calculation, especially when
     the system is not exactly at the equilibrium position. The 'accurate' choice
     is longer to execute (e.g. 3-6 times slower), but retains full forces. The
     'fast' mode also skips the IR and Raman mode calculation. The 'very fast'
     option halves the displacements, but assumes the initial configuration is
     fully equilibrated (which can be done with e.g. 'optimizer=BFGS'.
   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
     'PZ', 'PBE','revPBE','RPBE','PW91','VWN'     for Dacapo/Jacapo
     'LDA','PBE','revPBE','RPBE','GGA'            for ABINIT
     'LDA','PBE','REVPBE','PBESOL','WC06','AM05'  for ELK
     'LDA','PBE','RHF','MP2','B3LYP'              for NWChem
     '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-5 [eV].

 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' Type of database (ABINIT)
                 'sv','pv', ... (VASP)
   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)

 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.

 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).
 NWChem M. Valiev et al, Comput. Phys. Commun. 181, 1477 (2010).
   <http://www.nwchem-sw.org/>. Exists as Debian package 'nwchem' and 'nwchem-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
 DACAPO B. Hammer et al, Phys.Rev. B 59, 7413 (1999) 
   <https://wiki.fysik.dtu.dk/dacapo/>. Exists as Debian package 'dacapo' and 'dacapo-psp'.
 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.

 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: sqw_powder(model)
 To generate the dispersion curves along principal crystallographic directions,
   use: sqw_kpath(model)

 When performing a model evaluation, the DOS (phonon density of states) is also 
   computed and stored when the options 'dos' is specified during the creation. 
   The DOS is only computed during the first evaluation, and is stored in 
   model.UserData.DOS as an iData object. Subsequent evaluations are faster.

 In order to properly evaluate the neutron scattering intensity, 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 guess 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','dos');
   qh=linspace(0.01,.5,50);qk=qh; ql=qh; w=linspace(0.01,50,51);
   f=iData(s,[],qh,qk,ql',w); scatter3(log(f(1,:, :,:)),'filled');
   figure; plot(s.UserData.DOS); % plot the DOS, as indicated during model creation

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

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

 Version: Mar. 22, 2017
 See also iData, iFunc/fits, iFunc/plot, gauss, sqw_cubic_monoatomic, sqw_sine3d, sqw_vaks
   sqw_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 Wed 22-Mar-2017 09:13:30 by m2html © 2005. iFit (c) E.Farhi/ILL EUPL 1.1