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