1import sys 2import numpy as np 3 4from ase.units import _hbar, _c, _e, _me, Hartree, Bohr 5from gpaw import __version__ as version 6from gpaw.utilities.folder import Folder 7 8 9def get_folded_spectrum( 10 exlist=None, 11 emin=None, 12 emax=None, 13 de=None, 14 energyunit='eV', 15 folding='Gauss', 16 width=0.08, # Gauss/Lorentz width 17 form='r'): 18 """Return folded spectrum.""" 19 20 x = [] 21 y = [] 22 for ex in exlist: 23 x.append(ex.get_energy() * Hartree) 24 y.append(ex.get_oscillator_strength(form)) 25 26 if energyunit == 'nm': 27 # transform to experimentally used wavelength [nm] 28 x = 1.e+9 * 2 * np.pi * _hbar * _c / _e / np.array(x) 29 elif energyunit != 'eV': 30 raise RuntimeError('currently only eV and nm are supported') 31 32 return Folder(width, folding).fold(x, y, de, emin, emax) 33 34 35def spectrum(exlist=None, 36 filename=None, 37 emin=None, 38 emax=None, 39 de=None, 40 energyunit='eV', 41 folding='Gauss', 42 width=0.08, # Gauss/Lorentz width 43 comment=None, 44 form='r'): 45 """Write out a folded spectrum. 46 47 Parameters 48 ---------- 49 exlist: ExcitationList 50 filename: 51 File name for the output file, STDOUT if not given 52 emin: 53 min. energy, set to cover all energies if not given 54 emax: 55 max. energy, set to cover all energies if not given 56 de: 57 energy spacing 58 energyunit: {'eV', 'nm'} 59 Energy unit, default 'eV' 60 folding: {'Gauss', 'Lorentz'} 61 width: 62 folding width in terms of the chosen energyunit 63 """ 64 65 # output 66 out = sys.stdout 67 if filename is not None: 68 out = open(filename, 'w') 69 if comment: 70 print('#', comment, file=out) 71 72 print('# Photoabsorption spectrum from linear response TD-DFT', file=out) 73 print('# GPAW version:', version, file=out) 74 if folding is not None: # fold the spectrum 75 print('# %s folded, width=%g [%s]' % (folding, width, 76 energyunit), file=out) 77 if form == 'r': 78 out.write('# length form') 79 else: 80 assert(form == 'v') 81 out.write('# velocity form') 82 print('# om [%s] osz osz x osz y osz z' 83 % energyunit, file=out) 84 85 energies, values = get_folded_spectrum(exlist, emin, emax, de, 86 energyunit, folding, 87 width, form) 88 89 for e, val in zip(energies, values): 90 print('%10.5f %12.7e %12.7e %11.7e %11.7e' % 91 (e, val[0], val[1], val[2], val[3]), file=out) 92 93 if filename is not None: 94 out.close() 95 96 97def get_adsorbance_pre_factor(atoms): 98 """Return the absorbance pre-factor for solids. Unit m^-1. 99 100 Use this factor to multiply the folded oscillatior strength 101 obtained from spectrum() 102 103 Robert C. Hilborn, Am. J. Phys. 50, 982 (1982) 104 """ 105 return np.pi * _e**2 / 2. / _me / _c 106 107 108def rotatory_spectrum(exlist=None, 109 filename=None, 110 emin=None, 111 emax=None, 112 de=None, 113 energyunit='eV', 114 folding='Gauss', 115 width=0.08, # Gauss/Lorentz width 116 comment=None 117 ): 118 """Write out a folded rotatory spectrum. 119 120 See spectrum() for explanation of the parameters. 121 """ 122 123 # output 124 out = sys.stdout 125 if filename is not None: 126 out = open(filename, 'w') 127 if comment: 128 print('#', comment, file=out) 129 130 print('# Rotatory spectrum from linear response TD-DFT', file=out) 131 print('# GPAW version:', version, file=out) 132 if folding is not None: # fold the spectrum 133 print('# %s folded, width=%g [%s]' % (folding, width, 134 energyunit), file=out) 135 print('# om [%s] R [cgs]' 136 % energyunit, file=out) 137 138 x = [] 139 y = [] 140 for ex in exlist: 141 x.append(ex.get_energy() * Hartree) 142 y.append(ex.get_rotatory_strength()) 143 144 if energyunit == 'nm': 145 # transform to experimentally used wavelength [nm] 146 x = 1.e+9 * 2 * np.pi * _hbar * _c / _e / np.array(x) 147 y = np.array(y) 148 elif energyunit != 'eV': 149 raise RuntimeError('currently only eV and nm are supported') 150 151 energies, values = Folder(width, folding).fold(x, y, de, emin, emax) 152 for e, val in zip(energies, values): 153 print('%10.5f %12.7e' % (e, val), file=out) 154 155 if filename is not None: 156 out.close() 157 158 159class Writer(Folder): 160 """Object to write data to a file""" 161 def __init__(self, folding=None, width=0.08): 162 """Evaluate the polarizability from sum over states. 163 164 Parameters 165 ---------- 166 folding : string or None(default) 167 Name of the folding function, see gpaw/folder.py for options. 168 None means do not fold at all. 169 width : float 170 The width to be transferred to the folding function. 171 """ 172 self.folding = folding 173 if folding is not None: 174 Folder.__init__(self, width, folding) 175 176 def write(self, filename=None, 177 emin=None, emax=None, de=None, 178 comment=None): 179 180 out = sys.stdout 181 if filename is not None: 182 out = open(filename, 'w') 183 184 print('#', self.title, file=out) 185 print('# GPAW version:', version, file=out) 186 if comment: 187 print('#', comment, file=out) 188 if self.folding is not None: 189 print('# %s folded, width=%g [eV]' % (self.folding, 190 self.width), file=out) 191 energies, values = self.fold(self.energies, self.values, 192 de, emin, emax) 193 else: 194 energies, values = self.energies, self.values 195 196 print('#', self.fields, file=out) 197 198 for e, val in zip(energies, values): 199 string = '%10.5f' % e 200 for vf in val: 201 string += ' %12.7e' % vf 202 print(string, file=out) 203 204 if filename is not None: 205 out.close() 206 207 208def polarizability(exlist, omega, form='v', 209 tensor=False, index=0): 210 """Evaluate the polarizability from sum over states. 211 212 Parameters 213 ---------- 214 exlist: ExcitationList 215 omega: 216 Photon energy (eV) 217 form: {'v', 'r'} 218 Form of the dipole matrix element 219 index: {0, 1, 2, 3} 220 0: averaged, 1,2,3:alpha_xx, alpha_yy, alpha_zz 221 tensor: 222 if True returns alpha_ij, i,j=x,y,z 223 index is ignored 224 225 Returns 226 ------- 227 alpha: 228 Unit (e^2 Angstrom^2 / eV). 229 Multiply with Bohr * Ha to get (Angstrom^3) 230 """ 231 if tensor: 232 alpha = np.zeros(np.array(omega).shape + (3, 3), dtype=complex) 233 for ex in exlist: 234 alpha += ex.get_dipole_tensor(form=form) / ( 235 (ex.energy * Hartree)**2 - omega**2) 236 else: 237 alpha = np.zeros_like(omega) 238 for ex in exlist: 239 alpha += ex.get_oscillator_strength(form=form)[index] / ( 240 (ex.energy * Hartree)**2 - omega**2) 241 242 return alpha * Bohr**2 * Hartree 243