Branch data 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 : : !!***
|
