Line data Source code
1 : !!****m* ABINIT/m_model_screening
2 : !! NAME
3 : !! m_model_screening
4 : !!
5 : !! FUNCTION
6 : !! Module containing functions for calculating and fitting model dielectric functions
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MS)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_model_screening
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 :
28 : use m_io_tools, only : open_file
29 :
30 : implicit none
31 :
32 : private
33 :
34 : public :: im_screening ! Calc. Drude-Lorentz model function from parameters.
35 : public :: re_screening ! Calc. Drude-Lorentz model function from parameters.
36 : public :: re_and_im_screening ! Calc. Drude-Lorentz model function from parameters.
37 : public :: re_and_im_screening_with_phase ! Calc. Drude-Lorentz model function from parameters.
38 : ! with the addition of a phase
39 : public :: sequential_fitting ! Fit poles one by one
40 : public :: init_peaks_from_grid ! find approximate expression for parameters from
41 : ! chi0 or eps^-1 on a grid in the complex plane.
42 : public :: init_peaks_even_dist ! Initial guess from even distributuin of peaks
43 : public :: init_single_peak ! Initialise a single peak from the maximum, the
44 : ! origin, and the second value along the
45 : ! imaginary axis
46 : ! public :: int_screening ! Find the integral along real or complex axis
47 : ! from parameters.
48 : public :: remove_phase
49 :
50 : CONTAINS !==============================================================================
51 : !!***
52 :
53 : !!****f* m_model_screening/im_screening
54 : !! NAME
55 : !! im_screening
56 : !!
57 : !! FUNCTION
58 : !! Return the imaginary part of model dielectric / inverse dielectric
59 : !! function as given by a Drude-Lorentx model
60 : !!
61 : !! The function is a sum in the complex plane:
62 : !! f(z) = Sum_n f_n * Im[ ((w_n^2-z^2) - i*gamma*z)^-1 ], z=a-i*b
63 : !!
64 : !! Here, f_n is the oscillator strength, w_n the location of the peak for
65 : !! the imaginary function, and gamma is related to the width
66 : !!
67 : !! INPUTS
68 : !! omega = (complex) Real and imaginary part of the frequency points
69 : !! coeff = The coefficients in order: f_1,w_1,gamma_1,f_2,w_2,gamma_2,
70 : !! ...,f_n,w_n,gamma_n
71 : !! nomega = number of fit points
72 : !! ncoeff = number of coefficients
73 : !!
74 : !! OUTPUT
75 : !!
76 : !! NOTES
77 : !!
78 : !! SOURCE
79 :
80 0 : subroutine im_screening(omega,fval,nomega,coeff,ncoeff)
81 :
82 : !Arguments ------------------------------------
83 : !scalars
84 : integer,intent(in) :: nomega,ncoeff
85 : !arrays
86 : complex(dp),intent(in) :: omega(nomega)
87 : real(gwp) ,intent(in) :: coeff(ncoeff)
88 : real(gwp) ,intent(out) :: fval(nomega)
89 :
90 : !Local variables-------------------------------
91 : !scalars
92 : integer :: io,ip
93 : real(gwp) :: rez,imz,realp,imagp
94 : real(gwp) :: fn,wn,gamman
95 : ! *********************************************************************
96 :
97 : ! The expression is: -f_n*(2*rez*imz-rez*gamma_n)
98 : ! /( (-imz*gamma_n+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma_n)^2 )
99 :
100 0 : do io=1,nomega
101 0 : fval(io) = 0.0
102 0 : rez = REAL(omega(io))
103 0 : imz = AIMAG(omega(io))
104 0 : do ip=1,ncoeff,3
105 0 : fn = coeff(ip)
106 0 : wn = coeff(ip+1)
107 0 : gamman = coeff(ip+2)
108 0 : realp = -imz*gamman+wn*wn+imz*imz-rez*rez
109 0 : imagp = rez*(two*imz-gamman)
110 :
111 0 : fval(io) = fval(io)-fn*imagp/((realp*realp)+(imagp*imagp))
112 :
113 : end do
114 : end do
115 :
116 0 : end subroutine im_screening
117 : !!***
118 :
119 : !!****f* m_model_screening/re_screening
120 : !! NAME
121 : !! re_screening
122 : !!
123 : !! FUNCTION
124 : !! Return the real part of model dielectric / inverse dielectric
125 : !! function as evaluated from pole coefficients.
126 : !!
127 : !! The function is a sum of poles in the complex plane:
128 : !! f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
129 : !! where each pole occurs twice in a time-ordered fashion.
130 : !!
131 : !! Here, the A are the oscillator strengths, B the real component of the position
132 : !! of the pole, and C the imaginary component.
133 : !!
134 : !! INPUTS
135 : !! omega = (complex) Real and imaginary part of the frequency points
136 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
137 : !! nomega = number of fit points
138 : !! ncoeff = number of coefficients
139 : !!
140 : !! OUTPUT
141 : !!
142 : !! NOTES
143 : !!
144 : !! SOURCE
145 :
146 0 : subroutine re_screening(omega,fval,nomega,coeff,ncoeff)
147 :
148 : !Arguments ------------------------------------
149 : !scalars
150 : integer,intent(in) :: nomega,ncoeff
151 : !arrays
152 : complex(dp),intent(in) :: omega(nomega)
153 : real(gwp) ,intent(in) :: coeff(ncoeff)
154 : real(gwp) ,intent(out) :: fval(nomega)
155 :
156 : !Local variables-------------------------------
157 : !scalars
158 : integer :: io,ip
159 : real(gwp) :: rez,imz,realp,imagp
160 : real(gwp) :: fn,wn,gamman
161 : ! *********************************************************************
162 :
163 : ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
164 : ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
165 :
166 0 : do io=1,nomega
167 0 : fval(io) = 0.0
168 0 : rez = REAL(omega(io))
169 0 : imz = AIMAG(omega(io))
170 0 : do ip=1,ncoeff,3
171 0 : fn = coeff(ip)
172 0 : wn = coeff(ip+1)
173 0 : gamman = coeff(ip+2)
174 0 : realp = -imz*gamman+wn*wn+imz*imz-rez*rez
175 0 : imagp = rez*(two*imz-gamman)
176 :
177 0 : fval(io) = fval(io)-fn*realp/((realp*realp)+(imagp*imagp))
178 :
179 : end do
180 : end do
181 :
182 0 : end subroutine re_screening
183 : !!***
184 :
185 : !!****f* m_model_screening/re_and_im_screening
186 : !! NAME
187 : !! re_and_im_screening
188 : !!
189 : !! FUNCTION
190 : !! Return the real and imaginary part of model dielectric / inverse dielectric
191 : !! function as evaluated from pole coefficients.
192 : !!
193 : !! The function is a sum of poles in the complex plane:
194 : !! f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
195 : !! where each pole occurs twice in a time-ordered fashion.
196 : !!
197 : !! Here, the A are the oscillator strengths, B the real component of the position
198 : !! of the pole, and C the imaginary component.
199 : !!
200 : !! INPUTS
201 : !! omega = (complex) Real and imaginary part of the frequency points
202 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
203 : !! nomega = number of fit points
204 : !! ncoeff = number of coefficients
205 : !!
206 : !! OUTPUT
207 : !!
208 : !! NOTES
209 : !!
210 : !! SOURCE
211 :
212 0 : subroutine re_and_im_screening(omega,fval,nomega,coeff,ncoeff)
213 :
214 : !Arguments ------------------------------------
215 : !scalars
216 : integer,intent(in) :: nomega,ncoeff
217 : !arrays
218 : complex(dp) ,intent(in) :: omega(nomega)
219 : real(gwp) ,intent(in) :: coeff(ncoeff)
220 : complex(gwp),intent(out) :: fval(nomega)
221 :
222 : !Local variables-------------------------------
223 : !scalars
224 : integer :: io,ip
225 : real(gwp) :: rez,imz,realp,imagp,refval,imfval
226 : real(gwp) :: fn,wn,gamman
227 : ! *********************************************************************
228 :
229 : ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
230 : ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
231 :
232 0 : do io=1,nomega
233 0 : fval(io) = 0.0
234 0 : rez = REAL(omega(io))
235 0 : imz = AIMAG(omega(io))
236 0 : do ip=1,ncoeff,3
237 0 : fn = coeff(ip)
238 0 : wn = coeff(ip+1)
239 0 : gamman = coeff(ip+2)
240 0 : realp = -imz*gamman+wn*wn+imz*imz-rez*rez
241 0 : imagp = rez*(two*imz-gamman)
242 :
243 0 : refval = fn*realp/((realp*realp)+(imagp*imagp))
244 0 : imfval = fn*imagp/((realp*realp)+(imagp*imagp))
245 :
246 0 : fval(io) = fval(io)-CMPLX(refval,imfval)
247 :
248 : end do
249 : end do
250 :
251 0 : end subroutine re_and_im_screening
252 : !!***
253 :
254 : !!****f* m_model_screening/re_and_im_screening_with_phase
255 : !! NAME
256 : !! re_and_im_screening_with_phase
257 : !!
258 : !! FUNCTION
259 : !! Return the real and imaginary part of model dielectric / inverse dielectric
260 : !! function as evaluated from pole coefficients.
261 : !!
262 : !! The function is a sum of poles in the complex plane:
263 : !! f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
264 : !! where each pole occurs twice in a time-ordered fashion.
265 : !!
266 : !! Here, the A are the oscillator strengths, B the real component of the position
267 : !! of the pole, and C the imaginary component.
268 : !!
269 : !! INPUTS
270 : !! omega = (complex) Real and imaginary part of the frequency points
271 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
272 : !! nomega = number of fit points
273 : !! ncoeff = number of coefficients
274 : !!
275 : !! OUTPUT
276 : !!
277 : !! NOTES
278 : !!
279 : !! SOURCE
280 :
281 0 : subroutine re_and_im_screening_with_phase(omega,fval,nomega,coeff,ncoeff)
282 :
283 : !Arguments ------------------------------------
284 : !scalars
285 : integer,intent(in) :: nomega,ncoeff
286 : !arrays
287 : complex(dp) ,intent(in) :: omega(nomega)
288 : real(gwp) ,intent(in) :: coeff(ncoeff)
289 : complex(gwp),intent(out) :: fval(nomega)
290 :
291 : !Local variables-------------------------------
292 : !scalars
293 : integer :: io,ip,npoles
294 : real(gwp) :: rez,imz,realp,imagp,refval,imfval,retemp,imtemp
295 : real(gwp) :: fn,wn,gamman,imrot,rerot
296 : ! *********************************************************************
297 :
298 : ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
299 : ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
300 0 : npoles = (ncoeff-1)/3
301 :
302 0 : do io=1,nomega
303 0 : fval(io) = 0.0
304 0 : rez = REAL(omega(io))
305 0 : imz = AIMAG(omega(io))
306 0 : do ip=1,(ncoeff-1),3
307 0 : fn = coeff(ip)
308 0 : wn = coeff(ip+1)
309 0 : gamman = coeff(ip+2)
310 0 : realp = -imz*gamman+wn*wn+imz*imz-rez*rez
311 0 : imagp = rez*(two*imz-gamman)
312 :
313 0 : refval = fn*realp/((realp*realp)+(imagp*imagp))
314 0 : imfval = fn*imagp/((realp*realp)+(imagp*imagp))
315 :
316 0 : fval(io) = fval(io)-CMPLX(refval,imfval)
317 :
318 : end do
319 : ! Restore phase
320 0 : rerot = COS(coeff(npoles*3+1))
321 0 : imrot = SIN(coeff(npoles*3+1))
322 0 : retemp = REAL(fval(io))
323 0 : imtemp = AIMAG(fval(io))
324 0 : fval(io) = CMPLX(rerot*retemp-imrot*imtemp,rerot*imtemp + imrot*retemp)
325 : end do
326 :
327 0 : end subroutine re_and_im_screening_with_phase
328 : !!***
329 :
330 :
331 : !!****f* m_model_screening/sequential_fitting
332 : !! NAME
333 : !! sequential_fitting
334 : !!
335 : !! FUNCTION
336 : !! Fit a function in the complex plane pole-by-pole in such
337 : !! a way as to increasingly minimise the error
338 : !!
339 : !! INPUTS
340 : !! omega = (complex) Real and imaginary part of the frequency points
341 : !! refval = Real part of function to be fitted
342 : !! imfval = imaginary part of function to be fitted
343 : !! nomega = Total number of points in the complex plane
344 : !! nfreqre = Number of points along real axis
345 : !! nfreqim = Number of points along imaginary axis
346 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
347 : !! ncoeff = number of coefficients
348 : !! prtvol = Diagnostics verbose level
349 : !!
350 : !! OUTPUT
351 : !!
352 : !! NOTES
353 : !!
354 : !! SOURCE
355 :
356 0 : subroutine sequential_fitting(omega,refval,imfval,nomega,nfreqre,coeff,&
357 : & ncoeff,prtvol,startcoeff)
358 :
359 : !Arguments ------------------------------------
360 : !scalars
361 : integer,intent(in) :: nomega,nfreqre,ncoeff,prtvol
362 : !arrays
363 : complex(dp),intent(in) :: omega(nomega)
364 : real(gwp) ,intent(out) :: coeff(ncoeff)
365 : real(gwp) ,intent(inout) :: refval(nomega),imfval(nomega)
366 : real(gwp),optional,intent(out) :: startcoeff(ncoeff)
367 :
368 : !Local variables-------------------------------
369 : !scalars
370 : integer :: ip,npoles,idx
371 : real(gwp) :: thiscoeff(3),norm,invnorm
372 0 : real(dp) :: re_zvals(nomega),im_zvals(nomega)
373 : ! real(dp) :: orig_refval(nomega),orig_imfval(nomega)
374 0 : complex(gwp) :: pole_func(nomega)
375 : ! *********************************************************************
376 :
377 0 : npoles = ncoeff/3
378 0 : re_zvals(:) = REAL(omega(:))
379 0 : im_zvals(:) = AIMAG(omega(:))
380 :
381 : ! Normalise
382 0 : norm = MAXVAL(ABS(imfval))
383 0 : invnorm = 1.0_gwp/norm
384 0 : refval = invnorm*refval
385 0 : imfval = invnorm*imfval
386 :
387 : ! Loop over poles to fit
388 0 : do ip=1,npoles
389 0 : idx = 3*(ip-1)+1
390 : ! Initialise pole
391 0 : call init_single_peak(omega,refval,imfval,nomega,nfreqre,thiscoeff,prtvol)
392 0 : if (present(startcoeff)) then
393 0 : startcoeff(idx:idx+2) = thiscoeff(1:3)
394 : end if
395 : ! Make fit
396 : #ifdef HAVE_LEVMAR
397 : call dfit_re_and_im_screening(re_zvals,im_zvals,imfval,refval,&
398 : & nomega,3,thiscoeff,prtvol)
399 : #else
400 0 : ABI_ERROR(' ABINIT was not compiled with the levmar library!')
401 : #endif
402 : ! Remove current fit
403 0 : call re_and_im_screening(omega,pole_func,nomega,thiscoeff,3)
404 0 : refval(:) = refval(:) - REAL(pole_func(:))
405 0 : imfval(:) = imfval(:) - AIMAG(pole_func(:))
406 0 : coeff(idx:idx+2) = thiscoeff(1:3)
407 0 : coeff(idx) = norm*coeff(idx)
408 : end do
409 :
410 0 : end subroutine sequential_fitting
411 : !!***
412 :
413 : !!****f* m_model_screening/init_peaks_from_grid
414 : !! NAME
415 : !! init_peaks_from_grid
416 : !!
417 : !! FUNCTION
418 : !!
419 : !! Find an initial guess of coefficents from the "valleys" and "hills" in
420 : !! the complex plane.
421 : !!
422 : !! INPUTS
423 : !! omega = (complex) Real and imaginary part of the frequency points
424 : !! yvals = The function to be fitted in the complex plane
425 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
426 : !! nomega = number of fit points
427 : !! ncoeff = number of coefficients
428 : !! prtvol = Verbosity of diagnostics
429 : !!
430 : !! OUTPUT
431 : !!
432 : !! NOTES
433 : !!
434 : !! SOURCE
435 :
436 0 : subroutine init_peaks_from_grid(omega,fval,nomega,nfreqre,nfreqim,coeff,ncoeff,prtvol)
437 :
438 : !Arguments ------------------------------------
439 : !scalars
440 : integer,intent(in) :: nomega,nfreqre,nfreqim,ncoeff,prtvol
441 : !arrays
442 : complex(dp) ,intent(in) :: omega(nomega)
443 : real(gwp) ,intent(out) :: coeff(ncoeff)
444 : complex(gwp),intent(in) :: fval(nomega)
445 :
446 : !Local variables-------------------------------
447 : !scalars
448 : integer :: npoles,iline,idx,ip
449 : real(gwp) :: pol,gam,maxv,df,dk,val2,b2,osc,val1
450 : real(gwp) :: temp1,temp2,temp3
451 :
452 : !arrays
453 0 : integer :: ploc(ncoeff/3)
454 : ! *********************************************************************
455 :
456 0 : npoles = ncoeff/3
457 :
458 : ! Map the evolution of peaks if prtvol>10
459 0 : if (prtvol>10) then
460 0 : call print_peaks(omega,fval,nomega,nfreqre,nfreqim)
461 : end if ! prtvol>10
462 :
463 : ! Count the number of peaks per line and find the location of the
464 : ! constant-imaginary frequency line wich has at least a number of
465 : ! peaks commensurate with the requested number of poles
466 0 : call find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline)
467 0 : write(std_out,*) ' Optimum peak locations:',ploc
468 0 : write(std_out,*) ' on iline:',iline
469 : ! Now fit the peaks. A linear interpolation along the imaginary
470 : ! direction is used to get a rough estimate of the width of
471 : ! the peak.
472 0 : do ip=1,npoles
473 0 : pol = REAL(omega(ploc(ip)))
474 0 : maxv = AIMAG(fval(ploc(ip)))
475 0 : write(std_out,*) ' maxv:',maxv
476 0 : if (ploc(ip)<nfreqre+1) then ! We are right on the real axis
477 0 : if (ploc(ip)==1) then ! Peak is at origin (Drude peak, i.e. metal)
478 0 : b2 = AIMAG(omega(nfreqre+1))
479 0 : val2 = AIMAG(fval(nfreqre+1))
480 0 : write(std_out,*) '1: ploc:',ploc(ip),' b2:',b2,' val2:',val2
481 : else ! Second value will be in c-plane
482 0 : idx = nfreqre+nfreqim+ploc(ip)-1
483 0 : b2 = AIMAG(omega(idx))
484 0 : val2 = AIMAG(fval(idx))
485 0 : write(std_out,*) '2: ploc:',ploc(ip),' b2:',b2,' val2:',val2
486 : end if
487 0 : else if (ploc(ip)<nfreqre+nfreqim+1) then ! We are right on the imaginary axis
488 0 : if (ploc(ip)==nfreqre+nfreqim) then
489 0 : ABI_ERROR(' Peak in upper left corner. This should never happen')
490 : end if
491 0 : b2 = AIMAG(omega(ploc(ip)+1))
492 0 : val2 = AIMAG(fval(ploc(ip)+1))
493 0 : write(std_out,*) '3: ploc:',ploc(ip),' b2:',b2,' val2:',val2
494 : else ! We are in the complex plane
495 0 : idx = ploc(ip)+nfreqre-1
496 0 : b2 = AIMAG(omega(idx))
497 0 : val2 = AIMAG(fval(idx))
498 0 : write(std_out,*) '4: ploc:',ploc(ip),' idx:',idx,' b2:',b2,' val2:',val2
499 : end if
500 0 : df = ABS(val2 - maxv)
501 0 : dk = df/b2
502 0 : gam = -ABS(val2/dk)
503 : !temp1 = SQRT(-b2*b2*val2*(val2+maxv)+(maxv*maxv)**2)
504 : !temp2 = b2*(two*pol*pol+b2*b2)*val2-b2*maxv*pol*pol
505 : !temp3 = (pol*pol+b2*b2)*val2+maxv*pol*pol
506 : !gam = -((temp1-temp2)/temp3)
507 0 : if (gam>zero) gam = ((temp1+temp2)/temp3)
508 0 : osc = maxv*gam*pol
509 0 : idx = 3*(ip-1)+1
510 0 : coeff(idx ) = osc ! Oscillator strength
511 0 : coeff(idx+1) = pol ! Position of maximum
512 0 : coeff(idx+2) = gam ! Spread of function
513 0 : if (prtvol>9) then
514 0 : write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip
515 0 : write(std_out,'(a,ES16.8)') ' Osc. strength:',osc
516 0 : write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV'
517 0 : write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV'
518 0 : val2 = gam*half
519 0 : val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
520 0 : write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV'
521 0 : write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV'
522 0 : write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1)
523 : end if
524 : end do
525 :
526 0 : end subroutine init_peaks_from_grid
527 : !!***
528 :
529 : !!****f* m_model_screening/init_single_peak
530 : !! NAME
531 : !! init_single_peak
532 : !!
533 : !! FUNCTION
534 : !! Initialise a single peak by using the behaviour along the imaginary axis
535 : !! and the main peak
536 : !!
537 : !! INPUTS
538 : !! omega = (complex) Real and imaginary part of the frequency points
539 : !! yvals = The function to be fitted in the complex plane
540 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
541 : !! nomega = number of fit points
542 : !! ncoeff = number of coefficients
543 : !! prtvol = Verbosity of diagnostics
544 : !!
545 : !! OUTPUT
546 : !!
547 : !! NOTES
548 : !!
549 : !! SOURCE
550 :
551 0 : subroutine init_single_peak(omega,refval,imfval,nomega,nfreqre,coeff,prtvol)
552 :
553 : !Arguments ------------------------------------
554 : !scalars
555 : integer,intent(in) :: nomega,nfreqre,prtvol
556 : !arrays
557 : complex(dp),intent(in) :: omega(nomega)
558 : real(gwp) ,intent(out) :: coeff(3)
559 : real(gwp) ,intent(in) :: refval(nomega),imfval(nomega)
560 :
561 : !Local variables-------------------------------
562 : !scalars
563 : integer :: maxpos,idx
564 : real(gwp) :: pol,osc,gam,val1,val2,pol_sq
565 :
566 : ! *********************************************************************
567 :
568 0 : maxpos = MAXLOC(ABS(imfval(1:nfreqre)),1)
569 0 : if (maxpos==1) maxpos = MAXLOC(ABS(imfval(2:nfreqre)),1)
570 0 : pol = REAL(omega(maxpos))
571 0 : pol_sq = pol*pol
572 0 : osc = -refval(1)*pol_sq
573 0 : idx = nfreqre+1
574 0 : val2 = refval(idx)
575 0 : val1 = osc+val2*pol_sq+val2*AIMAG(omega(idx))*AIMAG(omega(idx))
576 0 : gam = -ABS(val1/(AIMAG(omega(idx))*val2))
577 :
578 0 : coeff(1) = osc
579 0 : coeff(2) = pol
580 0 : coeff(3) = gam
581 :
582 0 : if (prtvol>9) then
583 0 : write(std_out,'(a,ES16.8)') ' Osc. strength:',osc
584 0 : write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV'
585 0 : write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV'
586 0 : val2 = gam*half
587 0 : val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
588 0 : write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV'
589 0 : write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV'
590 0 : write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1)
591 : end if
592 :
593 0 : end subroutine init_single_peak
594 : !!***
595 :
596 : !!****f* m_model_screening/init_peaks_even_dist
597 : !! NAME
598 : !! init_peaks_even_dist
599 : !!
600 : !! FUNCTION
601 : !!
602 : !! Distribute the peaks evenly along a line in the complex plane and
603 : !! normalise.
604 : !!
605 : !! INPUTS
606 : !! omega = (complex) Real and imaginary part of the frequency points
607 : !! yvals = The function to be fitted in the complex plane
608 : !! coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n
609 : !! nomega = number of fit points
610 : !! ncoeff = number of coefficients
611 : !! prtvol = Verbosity of diagnostics
612 : !!
613 : !! OUTPUT
614 : !!
615 : !! NOTES
616 : !!
617 : !! SOURCE
618 :
619 0 : subroutine init_peaks_even_dist(omega,fval,nomega,nfreqre,coeff,ncoeff,prtvol)
620 :
621 : !Arguments ------------------------------------
622 : !scalars
623 : integer,intent(in) :: nomega,nfreqre,ncoeff,prtvol
624 : !arrays
625 : complex(dp) ,intent(in) :: omega(nomega)
626 : real(gwp) ,intent(out) :: coeff(ncoeff)
627 : complex(gwp),intent(in) :: fval(nomega)
628 :
629 : !Local variables-------------------------------
630 : !scalars
631 : integer :: npoles,ip,idx,iw
632 : real(gwp) :: delta,norm,div,val1,val2,osc,pol,gam
633 : ! *********************************************************************
634 :
635 0 : npoles = ncoeff/3
636 0 : div = real(npoles,gwp)
637 :
638 0 : delta = (omega(nfreqre)-omega(1))/(div+1.0_gwp)
639 : ! Integrate function along real axis (trapezoid rule) and have normalised
640 : ! oscillator strengths
641 0 : norm = fval(1)*half
642 0 : do iw=2,nfreqre-1
643 0 : norm = norm + fval(iw)
644 : end do
645 0 : norm = norm + fval(nfreqre)*half
646 0 : norm = norm*(omega(nfreqre)-omega(1))/real(nfreqre,gwp)
647 0 : norm = norm/div
648 :
649 0 : do ip=1,npoles
650 0 : idx = 3*(ip-1)+1
651 0 : pol = delta*ip ! Position of maximum
652 0 : gam = 0.1_gwp ! Spread of function
653 0 : val2 = gam*half
654 0 : val1 = SQRT(pol*pol-val2*val2)
655 0 : osc = norm*val1*two
656 0 : coeff(idx ) = osc!*(-1.0_gwp)**(ip-1)
657 0 : coeff(idx+1) = pol ! Position of maximum
658 0 : coeff(idx+2) = -gam ! Spread of function
659 0 : if (prtvol>9) then
660 0 : write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip
661 0 : write(std_out,'(a,ES16.8)') ' Osc. strength:',osc
662 0 : write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV'
663 0 : write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV'
664 0 : val2 = gam*half
665 0 : val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
666 0 : write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV'
667 0 : write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV'
668 0 : write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1)
669 : end if
670 : end do
671 :
672 0 : end subroutine init_peaks_even_dist
673 : !!***
674 :
675 : !!****f* m_model_screening/print_peaks
676 : !! NAME
677 : !! print_peaks
678 : !!
679 : !! FUNCTION
680 : !!
681 : !! Find and output the location of peaks on the grid in a file
682 : !!
683 : !! INPUTS
684 : !! omega = (complex) Real and imaginary part of the frequency points
685 : !! fvals = The function to be fitted in the complex plane
686 : !! nomega = number of fit points
687 : !! nfreqre = number or imaginary gridlines
688 : !! nfreqim = number of real gridlines
689 : !!
690 : !! OUTPUT
691 : !!
692 : !! NOTES
693 : !!
694 : !! SOURCE
695 :
696 0 : subroutine print_peaks(omega,fval,nomega,nfreqre,nfreqim)
697 :
698 : !Arguments ------------------------------------
699 : !scalars
700 : integer,intent(in) :: nomega,nfreqre,nfreqim
701 : !arrays
702 : complex(dp) ,intent(in) :: omega(nomega)
703 : complex(gwp),intent(in) :: fval(nomega)
704 :
705 : !Local variables-------------------------------
706 : !scalars
707 : integer :: ire,iim,unt_tmp
708 : integer :: idx1,idx2,idx3
709 : real(gwp) :: rez,imz,val1,val2,val3
710 : character(len=500) :: msg
711 : ! *********************************************************************
712 :
713 0 : if (open_file("grid_peak_tree.dat", msg, newunit=unt_tmp) /= 0) then
714 0 : ABI_ERROR(msg)
715 : end if
716 :
717 0 : do iim=nfreqim,1,-1
718 : ! write(std_out,*) ' iim:',iim
719 : ! Check first points
720 0 : idx1 = nfreqre+iim
721 0 : idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)
722 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
723 0 : val1 = AIMAG(fval(idx1))
724 0 : val2 = AIMAG(fval(idx2))
725 0 : if (ABS(val1)>ABS(val2)) then
726 0 : rez = REAL(omega(idx1))
727 0 : imz = AIMAG(omega(idx1))
728 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1
729 : end if
730 : ! Do all but the last
731 0 : do ire=1,nfreqre-4
732 0 : idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire
733 0 : idx2 = idx1+1
734 0 : idx3 = idx1+2
735 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3
736 0 : rez = REAL(omega(idx2))
737 0 : imz = AIMAG(omega(idx2))
738 0 : val1 = AIMAG(fval(idx1))
739 0 : val2 = AIMAG(fval(idx2))
740 0 : val3 = AIMAG(fval(idx3))
741 0 : if (((val1<val2).AND.(val2>val3))) then
742 0 : if (sign(1.0_gwp,val2)<zero) CYCLE
743 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
744 0 : else if (((val1>val2).AND.(val2<val3))) then
745 0 : if (sign(1.0_gwp,val2)>zero) CYCLE
746 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
747 : end if
748 : end do
749 : ! Check last point
750 0 : idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3
751 0 : idx2 = idx1 + 1
752 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
753 0 : rez = REAL(omega(idx2))
754 0 : imz = AIMAG(omega(idx2))
755 0 : val1 = AIMAG(fval(idx1))
756 0 : val2 = AIMAG(fval(idx2))
757 0 : if (ABS(val1)<ABS(val2)) then
758 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
759 : end if
760 : end do
761 : ! finally, do the purely real axis
762 : ! Check first points
763 0 : idx1 = 1; idx2 = 2
764 0 : val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2))
765 0 : if (ABS(val1)>ABS(val2)) then
766 0 : rez = REAL(omega(idx1))
767 0 : imz = AIMAG(omega(idx1))
768 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1
769 : end if
770 0 : do ire=2,nfreqre-3
771 0 : idx1 = ire; idx2 = idx1+1; idx3 = idx1+2
772 0 : rez = REAL(omega(idx2))
773 0 : imz = AIMAG(omega(idx2))
774 0 : val1 = AIMAG(fval(idx1))
775 0 : val2 = AIMAG(fval(idx2))
776 0 : val3 = AIMAG(fval(idx3))
777 0 : if (((val1<val2).AND.(val2>val3))) then
778 0 : if (sign(1.0_gwp,val2)<zero) CYCLE
779 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
780 0 : else if (((val1>val2).AND.(val2<val3))) then
781 0 : if (sign(1.0_gwp,val2)>zero) CYCLE
782 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
783 : end if
784 : end do
785 : ! Check last point
786 0 : idx1 = nfreqre-2
787 0 : idx2 = idx1 + 1
788 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
789 0 : rez = REAL(omega(idx2))
790 0 : imz = AIMAG(omega(idx2))
791 0 : val1 = AIMAG(fval(idx1))
792 0 : val2 = AIMAG(fval(idx2))
793 0 : if (ABS(val1)<ABS(val2)) then
794 0 : write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
795 : end if
796 :
797 0 : close(unt_tmp)
798 :
799 0 : end subroutine print_peaks
800 : !!***
801 :
802 : !!****f* m_model_screening/find_peaks
803 : !! NAME
804 : !! find_peaks
805 : !!
806 : !! FUNCTION
807 : !!
808 : !! Find the location of the highest peaks along gridlines starting at the real axis
809 : !! and then moving towards higher imaginary frequencies. Stop when enough
810 : !! peaks to satisfy the number of poles needed has been found
811 : !!
812 : !! INPUTS
813 : !! omega = (complex) Real and imaginary part of the frequency points
814 : !! fvals = The function to be fitted in the complex plane
815 : !! nomega = number of fit points
816 : !! nfreqre = number or imaginary gridlines
817 : !! nfreqim = number of real gridlines
818 : !!
819 : !! OUTPUT
820 : !!
821 : !! NOTES
822 : !!
823 : !! SOURCE
824 :
825 0 : subroutine find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline)
826 :
827 : !Arguments ------------------------------------
828 : !scalars
829 : integer,intent(in) :: nomega,nfreqre,nfreqim,npoles
830 : integer, intent(inout) :: iline
831 : !arrays
832 : integer ,intent(inout) :: ploc(npoles)
833 : complex(gwp), intent(in) :: fval(nomega)
834 :
835 : !Local variables-------------------------------
836 : !scalars
837 : integer :: ire,iim,ipoles
838 : integer :: idx1,idx2,idx3,ipol
839 : real(gwp) :: val1,val2,val3
840 : !arrays
841 0 : integer :: ploc_prev(npoles)
842 0 : real :: pval(npoles),pval_prev(npoles)
843 : ! *********************************************************************
844 :
845 0 : ploc=-1; ploc_prev=-1
846 0 : pval=zero; pval_prev=zero; ipol=1; ipoles=0
847 :
848 : ! First do a line along the real axis
849 0 : idx1 = 1; idx2 = 2
850 0 : val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2))
851 0 : if (ABS(val1)>ABS(val2)) then
852 0 : ipoles = ipoles + 1
853 0 : ploc(1)=idx1; pval(1)=val1
854 0 : write(std_out,*) ' pval:',pval
855 0 : write(std_out,*) ' ploc:',ploc
856 : end if
857 0 : do ire=2,nfreqre-3
858 0 : idx1 = ire; idx2 = idx1+1; idx3 = idx1+2
859 0 : val1 = AIMAG(fval(idx1))
860 0 : val2 = AIMAG(fval(idx2))
861 0 : val3 = AIMAG(fval(idx3))
862 0 : if (((val1<val2).AND.(val2>val3))) then
863 0 : if (sign(1.0_gwp,val2)<zero) CYCLE
864 0 : ipoles = ipoles + 1
865 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
866 0 : ipol = MINLOC(ABS(pval(:)),1)
867 0 : ploc(ipol)=idx2; pval(ipol)=val2
868 0 : write(std_out,*) ' pval:',pval
869 0 : write(std_out,*) ' ploc:',ploc
870 : end if
871 0 : else if (((val1>val2).AND.(val2<val3))) then
872 0 : if (sign(1.0_gwp,val2)>zero) CYCLE
873 0 : ipoles = ipoles + 1
874 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
875 0 : ipol = MINLOC(ABS(pval(:)),1)
876 0 : ploc(ipol)=idx2; pval(ipol)=val2
877 0 : write(std_out,*) ' pval:',pval
878 0 : write(std_out,*) ' ploc:',ploc
879 : end if
880 : end if
881 : end do
882 : ! Check last point
883 0 : idx1 = nfreqre-2
884 0 : idx2 = idx1 + 1
885 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
886 0 : val1 = AIMAG(fval(idx1))
887 0 : val2 = AIMAG(fval(idx2))
888 0 : if (ABS(val1)<ABS(val2)) then
889 0 : ipoles = ipoles + 1
890 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
891 0 : ipol = MINLOC(ABS(pval(:)),1)
892 0 : ploc(ipol)=idx2; pval(ipol)=val2
893 0 : write(std_out,*) ' pval:',pval
894 0 : write(std_out,*) ' ploc:',ploc
895 : end if
896 : end if
897 0 : write(std_out,'(a,i0)') ' Number of poles real axis:',ipoles
898 :
899 :
900 0 : ploc_prev = ploc; pval_prev = pval
901 :
902 : ! Do the rest of the imaginary grid until total
903 : ! number of peaks found equals npoles or less
904 0 : do iim=1,nfreqim-1
905 0 : ploc=-1; pval=zero; ipol=1; ipoles=0
906 0 : idx1 = nfreqre+iim
907 0 : idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)
908 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
909 0 : val1 = AIMAG(fval(idx1))
910 0 : val2 = AIMAG(fval(idx2))
911 0 : if (ABS(val1)>ABS(val2)) then
912 0 : ipoles = ipoles + 1
913 0 : ploc(1)=idx1; pval(1)=val1
914 0 : write(std_out,*) ' pval:',pval
915 0 : write(std_out,*) ' ploc:',ploc
916 : end if
917 : ! Do all but the last
918 0 : do ire=1,nfreqre-4
919 0 : idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire
920 0 : idx2 = idx1+1
921 0 : idx3 = idx1+2
922 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3
923 0 : val1 = AIMAG(fval(idx1))
924 0 : val2 = AIMAG(fval(idx2))
925 0 : val3 = AIMAG(fval(idx3))
926 0 : if (((val1<val2).AND.(val2>val3))) then
927 0 : if (sign(1.0_gwp,val2)<zero) CYCLE
928 0 : ipoles = ipoles + 1
929 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
930 0 : ipol = MINLOC(ABS(pval(:)),1)
931 0 : ploc(ipol)=idx2; pval(ipol)=val2
932 0 : write(std_out,*) ' pval:',pval
933 0 : write(std_out,*) ' ploc:',ploc
934 : end if
935 0 : else if (((val1>val2).AND.(val2<val3))) then
936 0 : if (sign(1.0_gwp,val2)>zero) CYCLE
937 0 : ipoles = ipoles + 1
938 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
939 0 : ipol = MINLOC(ABS(pval(:)),1)
940 0 : ploc(ipol)=idx2; pval(ipol)=val2
941 0 : write(std_out,*) ' pval:',pval
942 0 : write(std_out,*) ' ploc:',ploc
943 : end if
944 : end if
945 : end do
946 : ! Check last point
947 0 : idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3
948 0 : idx2 = idx1 + 1
949 : ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2
950 0 : val1 = AIMAG(fval(idx1))
951 0 : val2 = AIMAG(fval(idx2))
952 0 : if (ABS(val1)<ABS(val2)) then
953 0 : ipoles = ipoles + 1
954 0 : if (ANY( ABS(pval(:))<ABS(val2) )) then
955 0 : ipol = MINLOC(ABS(pval(:)),1)
956 0 : ploc(ipol)=idx2; pval(ipol)=val2
957 0 : write(std_out,*) ' pval:',pval
958 0 : write(std_out,*) ' ploc:',ploc
959 : end if
960 : end if
961 0 : write(std_out,'(2(a,i0))') ' Line,:',iim,' ipoles:',ipoles
962 0 : if (ipoles<=npoles) then
963 0 : iline = iim - 1
964 0 : ploc = ploc_prev
965 : EXIT
966 : end if
967 :
968 0 : ploc_prev = ploc; pval_prev = pval
969 :
970 : end do
971 :
972 0 : end subroutine find_peaks
973 : !!***
974 :
975 : !!****f* m_model_screening/remove_phase
976 : !! NAME
977 : !! remove_phase
978 : !!
979 : !! FUNCTION
980 : !! Find out what the complex phase factor is for off-diagonal elements
981 : !! and unmix the components.
982 : !!
983 : !! INPUTS
984 : !! fvals = The function to be fitted in the complex plane.
985 : !! nomega = Number of fit points.
986 : !! nfreqre = number or imaginary gridlines
987 : !! phase = The phase angle.
988 : !!
989 : !! OUTPUT
990 : !!
991 : !! NOTES
992 : !!
993 : !! SOURCE
994 :
995 0 : subroutine remove_phase(fval,nomega,phase)
996 :
997 : !Arguments ------------------------------------
998 : !scalars
999 : integer, intent(in) :: nomega
1000 : real(gwp), intent(out) :: phase
1001 : !arrays
1002 : complex(gwp), intent(inout) :: fval(nomega)
1003 :
1004 : !Local variables-------------------------------
1005 : !scalars
1006 : integer :: io
1007 : real(gwp) :: a,b,retemp,imtemp
1008 : ! *********************************************************************
1009 :
1010 : ! The phase can be found by checking when the function is
1011 : ! identically zero along the imaginary axis
1012 0 : if (ABS(AIMAG(fval(1)))<tol14) then ! Phase is zero
1013 0 : phase = zero
1014 0 : RETURN
1015 0 : else if (ABS(REAL(fval(1)))<tol14) then ! Phase is exactly pi/2
1016 0 : phase = pi*half
1017 0 : a = zero
1018 0 : b = -1.0_gwp
1019 : else
1020 0 : phase = ATAN(AIMAG(fval(1))/REAL(fval(1)))
1021 0 : a = COS(phase)
1022 0 : b = -SIN(phase)
1023 : end if
1024 : ! Rotate values
1025 0 : do io=1,nomega
1026 0 : retemp = REAL(fval(io))
1027 0 : imtemp = AIMAG(fval(io))
1028 0 : fval(io) = CMPLX(a*retemp-b*imtemp,a*imtemp+b*retemp)
1029 : end do
1030 :
1031 : end subroutine remove_phase
1032 : !!***
1033 :
1034 : END MODULE m_model_screening
1035 : !!***
|
