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 OCTOPUS QuantumEspresso Plane-wave pseudopotential code SIESTA 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') sqw_phonons('gui') 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)' 'molecule("H2O")' '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> <https://wiki.fysik.dtu.dk/ase/tutorials/spacegroup/spacegroup.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: options.python='python3' 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: <http://phonondb.mtl.kyoto-u.ac.jp/> <https://www.materialsproject.org/> <http://nomad-repository.eu/cms/> WARNING: 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). <http://chianti.geol.ucl.ac.uk/~dario>. Quantum Espresso P. Giannozzi, et al, J.Phys.:Condens.Matter, 21, 395502 (2009). <http://www.quantum-espresso.org/>. 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) <https://atztogo.github.io/phonopy/> 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 intensity. 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(1)=Amplitude 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.