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)' '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> '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: <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. 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). <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. 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(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','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.