Line data Source code
1 : !!****m* ABINIT/m_paw_gaussfit
2 : !! NAME
3 : !! m_paw_gaussfit
4 : !!
5 : !! FUNCTION
6 : !! Module to fit PAW related data to sums of gaussians
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2012-2026 ABINIT group (T. Rangel)
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 : !! NOTES
15 : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
16 : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
17 : !!
18 : !! SOURCE
19 :
20 : #include "libpaw.h"
21 :
22 : module m_paw_gaussfit
23 :
24 : USE_DEFS
25 : USE_MSG_HANDLING
26 : USE_MPI_WRAPPERS
27 : USE_MEMORY_PROFILING
28 :
29 : use m_paw_numeric, only : paw_splint, paw_spline
30 : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_deducer0, pawrad_free, pawrad_ifromr
31 :
32 : implicit none
33 :
34 : private
35 :
36 : public:: gaussfit_projector !fit non-local projectors to a sum of gaussians
37 : public:: gaussfit_main !main routine to fit
38 : !Routines related to MPI:
39 : private:: gaussfit_mpi_set_weight
40 : private:: gaussfit_mpi_remove_item
41 : private:: gaussfit_mpi_add_item
42 : private:: gaussfit_mpi_assign
43 : private:: gaussfit_mpi_main
44 : private:: gaussfit_mpi_calc_deviation
45 : private:: gaussfit_mpi_swap
46 : !Routines related to fitting:
47 : private:: gaussfit_fit !fit a function to gaussians
48 : private:: gaussfit_calc_deriv_r !calc. derivatives for real gauss.
49 : private:: gaussfit_calc_deriv_c !calc. derivatives for cplex. gauss. of form 1
50 : private:: gaussfit_calc_deriv_c2 !calc. derivatives for cplex. gauss. of form 2
51 : private:: gaussfit_calc_deriv_c3 !calc. derivatives for cplex. gauss. of form 3
52 : private:: gaussfit_calc_deriv_c4 !calc. derivatives for cplex. gauss. of form 4
53 : private:: gaussfit_rlsf !retreined least squares fit
54 : private:: gaussfit_chisq_alpha_beta
55 : !set parameters for LSF (for 5 different forms of gauss. sums):
56 : private:: gaussfit_set_param1
57 : private:: gaussfit_set_param2
58 : private:: gaussfit_set_param3
59 : private:: gaussfit_set_param4
60 : private:: gaussfit_set_param5
61 : private:: gaussfit_constrains_init !initialize constrains
62 : private:: gaussfit_apply_constrains !apply cons. to get new params.
63 :
64 : !!***
65 : integer,private,parameter:: positive=2
66 : integer,private,parameter:: restricted=3
67 : integer,private,parameter:: restricted_and_positive=4
68 :
69 : CONTAINS
70 : !===========================================================
71 : !!***
72 :
73 : !!****f* m_paw_gaussfit/gaussfit_main
74 : !! NAME
75 : !! gaussfit_main
76 : !!
77 : !! FUNCTION
78 : !! Fits a given input function f(r) to a sum of gaussians
79 : !!
80 : !! INPUTS
81 : !! nterm_bounds= sets the minimum and maximum number of terms to be taken into account.
82 : !! mparam= maximum number of parameters
83 : !! if(option==1)mparam=nterm_bounds(2)*4
84 : !! if(option==2)mparam=nterm_bounds(2)*6
85 : !! if(option==3)mparam=nterm_bounds(2)*2
86 : !! if(option==4)mparam=nterm_bounds(2)*4
87 : !! nr= number of real space points
88 : !! pawrad= pawrad type
89 : !! option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
90 : !! 2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
91 : !! 3 fit to a1 cos (k x^2) + a2 sin (k x^2)
92 : !! 4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
93 : !! Given these definitions, the number of complex gaussians are:
94 : !! ngauss=4*nterm (for option=1,2)
95 : !! ngauss=2*nterm (for option=3,4)
96 : !! outfile= filename to write out fitted functions.
97 : !! rpaw=paw radius
98 : !! y(nr)= function to fit
99 : !!
100 : !! OUTPUT
101 : !! nparam_out= number of parameters found.
102 : !! param_out(nparam_out)= parameters (coefficients and factors of complex gaussians).
103 : !!
104 : !! SOURCE
105 :
106 0 : subroutine gaussfit_main(mparam,nparam_out,nterm_bounds,nr,&
107 0 : & param_out,pawrad,option,outfile,rpaw,y,comm_mpi)
108 :
109 : !Arguments ------------------------------------
110 : integer,intent(in)::mparam,nr,nterm_bounds(2)
111 : integer,intent(in)::option
112 : integer,intent(in),optional:: comm_mpi
113 : real(dp),intent(in)::rpaw
114 : real(dp),intent(inout)::y(nr)
115 : character(80),intent(in)::outfile
116 : type(pawrad_type),intent(in) :: pawrad
117 : integer,intent(out)::nparam_out
118 : real(dp),intent(out)::param_out(mparam)
119 :
120 : !Local variables-------------------------------
121 : !scalars
122 : logical,parameter::modify_y=.false. !used only for plotting purposes
123 : integer :: ichisq,ierr,ii,jj
124 : integer :: master,me
125 : integer :: maxiter,minterm,my_chisq_size,counts_all
126 : integer :: ngauss,nparam,nproc,nterm
127 : integer :: verbosity
128 : real(dp) :: chisq,chisq_min
129 : !real(dp) :: T1,T2 !uncomment for timming
130 : !arrays
131 0 : integer::constrains(mparam)
132 0 : integer::proc_dist(nterm_bounds(1):nterm_bounds(2))
133 0 : integer,allocatable::counts(:),disp(:),map_nterm(:)
134 0 : real(dp)::limit(mparam),param_tmp(mparam),weight(mparam)
135 0 : real(dp),allocatable::chisq_array(:),recv_buf(:),send_buf(:)
136 0 : real(dp),allocatable::y_out(:)
137 : character(len=500) :: msg
138 :
139 : ! *************************************************************************
140 :
141 : !initialize variables
142 0 : maxiter=200
143 : !
144 : !initialize mpi quantities:
145 0 : master=0; me=0; nproc=1; proc_dist=1;
146 0 : if(present(comm_mpi)) then
147 0 : me=xmpi_comm_rank(comm_mpi)
148 0 : nproc=xmpi_comm_size(comm_mpi)
149 : end if
150 0 : if(nproc>1) then
151 : ! Find distribution (master)
152 0 : if(me==master) then
153 0 : call gaussfit_mpi_main(nproc,nterm_bounds,proc_dist)
154 : end if
155 : ! send distribution to all processors
156 : call xmpi_bcast(proc_dist(nterm_bounds(1):nterm_bounds(2)),&
157 0 : & master,comm_mpi,ierr)
158 : end if
159 : !Set size of chisq treated by each proc
160 0 : my_chisq_size=nterm_bounds(2)-nterm_bounds(1)+1 !all terms
161 0 : if(nproc>1 .and. .not. me==master) then
162 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
163 0 : if(.not. proc_dist(nterm)==me+1) cycle
164 0 : my_chisq_size=my_chisq_size+1
165 : end do
166 : end if
167 :
168 : !
169 : !Allocate objects
170 : !
171 0 : LIBPAW_ALLOCATE(y_out,(nr))
172 0 : LIBPAW_ALLOCATE(chisq_array,(my_chisq_size))
173 0 : if(master==me ) then
174 0 : LIBPAW_BOUND1_ALLOCATE(map_nterm,BOUNDS(nterm_bounds(1),nterm_bounds(2)))
175 0 : jj=1
176 0 : do ii=1,nproc
177 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
178 0 : if(proc_dist(nterm)==ii) then
179 0 : map_nterm(nterm)=jj
180 0 : jj=jj+1
181 : end if
182 : end do
183 : end do
184 : end if
185 : !
186 : !fill with zeros
187 : !
188 0 : chisq_array=zero
189 0 : nparam_out=0
190 0 : y_out=0.d0
191 0 : verbosity=1 !print the minimum at the terminal
192 : !
193 0 : ichisq=0
194 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
195 : ! mpi distribution
196 0 : if(.not. proc_dist(nterm)==me+1) cycle
197 0 : ichisq=ichisq+1
198 :
199 : ! call CPU_TIME(T1)
200 :
201 0 : if(option==1) nparam=nterm*4
202 0 : if(option==2) nparam=nterm*6
203 0 : if(option==3) nparam=nterm*2
204 0 : if(option==4) nparam=nterm*4
205 : ! set initial guess
206 : ! if(option==1) then
207 : ! call gaussfit_set_param3(nterm,nparam,param_tmp(1:nparam),sep(ii))
208 : ! elseif(option==2) then
209 : ! call gaussfit_set_param1(nterm,nparam,nr,&
210 : !& param_tmp(1:nparam),sep(ii),pawrad%rad(1:nr),y)
211 0 : if(option==3) then
212 0 : call gaussfit_set_param4(nparam,param_tmp(1:nparam))
213 0 : elseif(option==4) then
214 : call gaussfit_set_param5(nterm,nparam,nr,&
215 0 : & param_tmp(1:nparam),rpaw,y)
216 : end if
217 : !
218 : !
219 : call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),&
220 0 : & limit(1:nparam),nparam,nterm,nr,option,rpaw,y)
221 : !
222 : call gaussfit_fit(chisq,constrains(1:nparam),&
223 : & limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,&
224 0 : & verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out)
225 :
226 : ! if there was an error, set chisq to very high
227 : ! if there is an NaN, for instance:
228 0 : if(abs(chisq+1.d0)<tol8) chisq=99999
229 0 : if(abs(chisq)==chisq*chisq) chisq=99999
230 0 : if(chisq .ne. chisq) chisq=99999
231 : !
232 0 : chisq_array(ichisq)=chisq
233 :
234 : ! call CPU_TIME(T2)
235 : ! print *, 'Time for fit ', T2-T1, 'seconds.'
236 :
237 : end do !nterm
238 :
239 : !mpicast results
240 : !send distribution to master
241 0 : if(nproc>1) then
242 : ! Prepare communications:
243 0 : LIBPAW_ALLOCATE(counts,(nproc))
244 0 : counts(:)=0
245 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
246 0 : counts(proc_dist(nterm))=counts(proc_dist(nterm))+1
247 : end do
248 0 : counts_all=sum(counts)
249 0 : LIBPAW_ALLOCATE(send_buf,(counts(me+1)))
250 0 : send_buf(:)=chisq_array(1:counts(me+1))
251 0 : if(me==master) then
252 0 : LIBPAW_ALLOCATE(recv_buf,(counts_all))
253 : else
254 0 : LIBPAW_ALLOCATE(recv_buf,(1))
255 : end if
256 0 : LIBPAW_ALLOCATE(disp,(nproc))
257 0 : disp(1)=0
258 0 : do ii=2,nproc
259 0 : disp(ii)=disp(ii-1)+counts(ii-1)
260 : end do
261 : ! communicate all info to master
262 : call xmpi_gatherv(send_buf,counts(me+1),recv_buf,&
263 0 : & counts,disp,master,comm_mpi,ierr)
264 : ! fill in chisq_array with all received info:
265 0 : if(master==me) then
266 0 : do ii=1,counts_all
267 0 : chisq_array(ii)=recv_buf(ii)
268 : end do
269 : end if
270 : ! Deallocate MPI arrays:
271 0 : LIBPAW_DEALLOCATE(recv_buf)
272 0 : LIBPAW_DEALLOCATE(counts)
273 0 : LIBPAW_DEALLOCATE(disp)
274 0 : LIBPAW_DEALLOCATE(send_buf)
275 : end if
276 :
277 : !Print out info:
278 0 : if(me==master) then
279 0 : write(msg,'(3a)')'Preliminary results (with only 200 iter.):',ch10,' ngauss chisq'
280 0 : call wrtout(std_out,msg,'COLL')
281 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
282 0 : if(option==1) ngauss=nterm*4
283 0 : if(option==2) ngauss=nterm*4
284 0 : if(option==3) ngauss=nterm*2
285 0 : if(option==4) ngauss=nterm*2
286 0 : write(msg,'(i4,2x,e13.6,1x)')ngauss,&
287 0 : & chisq_array(map_nterm(nterm))
288 0 : call wrtout(std_out,msg,'COLL')
289 : end do
290 : end if
291 :
292 : !get minterm for best accuracy:
293 0 : if(me==master) then
294 0 : chisq_min=999999999
295 0 : do nterm=nterm_bounds(1),nterm_bounds(2)
296 0 : if(chisq_array(map_nterm(nterm))<chisq_min) then
297 0 : chisq_min=chisq_array(map_nterm(nterm))
298 0 : minterm=nterm
299 : end if
300 : end do
301 :
302 : !run again with the best parameters
303 0 : nterm=minterm
304 0 : if(option==1)nparam=4*nterm
305 0 : if(option==2)nparam=6*nterm
306 0 : if(option==3)nparam=2*nterm
307 0 : if(option==4)nparam=4*nterm
308 :
309 : ! set initial guess
310 0 : if (option==3) then
311 0 : call gaussfit_set_param4(nparam,param_tmp(1:nparam))
312 0 : elseif (option==4) then
313 : call gaussfit_set_param5(nterm,nparam,nr,&
314 0 : & param_tmp(1:nparam),rpaw,y)
315 : end if
316 : !
317 : call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),&
318 0 : & limit(1:nparam),nparam,nterm,nr,option,rpaw,y)
319 : !
320 : verbosity=1;
321 0 : maxiter=1000 !this time we do more iterations
322 : call gaussfit_fit(chisq,constrains(1:nparam),&
323 : & limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,&
324 0 : & verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out)
325 :
326 : ! Write out best solution
327 0 : write(msg,'(3a)')"Best solution (with more iterations):",ch10," ngauss chisq"
328 0 : call wrtout(std_out,msg,'COLL')
329 0 : if(option==1) ngauss=nterm*4
330 0 : if(option==2) ngauss=nterm*4
331 0 : if(option==3) ngauss=nterm*2
332 0 : if(option==4) ngauss=nterm*2
333 0 : write(msg,'(i4,2x,e13.6,1x)')ngauss,chisq
334 0 : call wrtout(std_out,msg,'COLL')
335 : end if
336 :
337 : !Fill output variables
338 : !and communicate results to all procs:
339 0 : if(option==1) then
340 0 : nparam_out=minterm*4
341 0 : elseif (option==2) then
342 0 : nparam_out=minterm*6
343 0 : elseif (option==3) then
344 0 : nparam_out=minterm*2
345 0 : elseif (option==4) then
346 0 : nparam_out=minterm*4
347 : end if
348 :
349 : !communicate
350 0 : if(nproc>1) then
351 0 : call xmpi_bcast(nparam_out,master,comm_mpi,ierr)
352 0 : call xmpi_bcast(param_tmp(1:nparam_out),master,comm_mpi,ierr)
353 : end if !nproc>1
354 :
355 0 : param_out(:)=param_tmp
356 :
357 : if(modify_y) then
358 : ! call xmpi_scatterv(y_out,nr,mpi_displ,y_out,nr,0,comm_mpi,ierr)
359 : y=y_out !at output modify y for the fitted y
360 : end if
361 :
362 0 : LIBPAW_DEALLOCATE(y_out)
363 0 : LIBPAW_DEALLOCATE(chisq_array)
364 0 : if(me==master) then
365 0 : LIBPAW_DEALLOCATE(map_nterm)
366 : end if
367 :
368 0 : end subroutine gaussfit_main
369 : !!***
370 :
371 : !----------------------------------------------------------------------
372 :
373 : !!****f* m_paw_gaussfit/gaussfit_mpi_set_weight
374 : !! NAME
375 : !! gaussfit_mpi_set_weight
376 : !!
377 : !! FUNCTION
378 : !! It sets a weight to the number of
379 : !! Gaussians used.
380 : !! This was calculated by measuring the time
381 : !! it takes to fit a projector with different
382 : !! number of gaussians.
383 : !!
384 : !! INPUTS
385 : !!
386 : !! OUTPUT
387 : !!
388 : !! SOURCE
389 :
390 0 : subroutine gaussfit_mpi_set_weight(f,x)
391 :
392 : !Arguments ------------------------------------
393 : integer,intent(in)::x
394 : integer,intent(out)::f
395 :
396 : !Local variables ------------------------------
397 : real(dp)::a,b,c,d,ff,xx
398 :
399 : !************************************************************************
400 :
401 : !The following parameters were obtained
402 : !from the time (in seconds)
403 : !it takes to fit a given projector
404 : !using from 1 to 100 gaussians.
405 0 : a = -0.374137d0 ! +/- 2.013 (538.1%)
406 0 : b = 0.207854d0 ! +/- 0.3385 (162.8%)
407 0 : c = 0.0266371d0 ! +/- 0.01534 (57.59%)
408 0 : d = 0.000152476d0 ! +/- 0.0001978 (129.7%)
409 :
410 0 : xx=real(x,dp)
411 0 : ff=a+b*xx+c*xx**2+d*xx**3
412 0 : f=max(1,ceiling(ff))
413 :
414 0 : end subroutine gaussfit_mpi_set_weight
415 : !!***
416 :
417 : !----------------------------------------------------------------------
418 :
419 : !!****f* m_paw_gaussfit/gaussfit_mpi_remove_item
420 : !! NAME
421 : !! gaussfit_mpi_remove_item
422 : !!
423 : !! FUNCTION
424 : !!
425 : !! INPUTS
426 : !!
427 : !! OUTPUT
428 : !!
429 : !! SOURCE
430 :
431 0 : subroutine gaussfit_mpi_remove_item(iterm,pload)
432 :
433 : !Arguments ------------------------------------
434 : integer,intent(in)::iterm
435 : integer,intent(inout)::pload
436 :
437 : !Local variables ------------------------------
438 : integer:: f_i
439 :
440 : !***********************************************************************
441 :
442 0 : call gaussfit_mpi_set_weight(f_i,iterm)
443 0 : pload=pload-f_i
444 :
445 0 : end subroutine gaussfit_mpi_remove_item
446 : !!***
447 :
448 : !----------------------------------------------------------------------
449 :
450 : !!****f* m_paw_gaussfit/gaussfit_mpi_add_item
451 : !! NAME
452 : !! gaussfit_mpi_add_item
453 : !!
454 : !! FUNCTION
455 : !!
456 : !! INPUTS
457 : !!
458 : !! OUTPUT
459 : !!
460 : !! SOURCE
461 :
462 0 : subroutine gaussfit_mpi_add_item(iterm,pload)
463 :
464 : !Arguments ------------------------------------
465 : integer,intent(in)::iterm
466 : integer,intent(inout)::pload
467 :
468 : !Local variables ------------------------------
469 : integer:: f_i
470 :
471 : !************************************************************************
472 :
473 0 : call gaussfit_mpi_set_weight(f_i,iterm)
474 0 : pload=pload+f_i
475 :
476 : end subroutine gaussfit_mpi_add_item
477 : !!***
478 :
479 : !----------------------------------------------------------------------
480 :
481 : !!****f* m_paw_gaussfit/gaussfit_mpi_calc_deviation
482 : !! NAME
483 : !! gaussfit_mpi_calc_deviation
484 : !!
485 : !! FUNCTION
486 : !!
487 : !! INPUTS
488 : !!
489 : !! OUTPUT
490 : !!
491 : !! SOURCE
492 :
493 0 : subroutine gaussfit_mpi_calc_deviation(deviation,nproc,proc_load)
494 :
495 : !Arguments ------------------------------------
496 : integer,intent(in)::nproc
497 : integer,intent(in)::proc_load(nproc)
498 : integer,intent(out)::deviation
499 :
500 : !Local variables ------------------------------
501 : integer:: jproc,kproc,kload,jload
502 :
503 : !************************************************************************
504 :
505 : ! deviation=0
506 : ! do jproc=1,nproc
507 : ! deviation=deviation+abs(proc_load(jproc)-ideal)
508 : ! end do
509 :
510 0 : deviation=0
511 0 : do jproc=1,nproc
512 0 : jload=proc_load(jproc)
513 0 : do kproc=1,jproc
514 0 : kload=proc_load(kproc)
515 0 : deviation=deviation+abs(kload-jload)
516 : end do
517 : end do
518 :
519 0 : end subroutine gaussfit_mpi_calc_deviation
520 : !!***
521 :
522 : !----------------------------------------------------------------------
523 :
524 : !!****f* m_paw_gaussfit/gaussfit_mpi_swap
525 : !! NAME
526 : !! gaussfit_mpi_swap
527 : !!
528 : !! FUNCTION
529 : !!
530 : !! INPUTS
531 : !!
532 : !! OUTPUT
533 : !!
534 : !! SOURCE
535 :
536 0 : subroutine gaussfit_mpi_swap(iterm,jterm,&
537 0 : & nproc,nterm_bounds,proc_dist,proc_load)
538 :
539 : !Arguments ------------------------------------
540 : integer,intent(in)::iterm,jterm,nproc,nterm_bounds(2)
541 : integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc)
542 :
543 : !Local variables ------------------------------
544 : integer:: deviation1,deviation2
545 : integer:: iproc,jproc
546 :
547 : !************************************************************************
548 :
549 : !Calculate initial state
550 0 : call gaussfit_mpi_calc_deviation(deviation1,nproc,proc_load)
551 0 : iproc=proc_dist(iterm)
552 0 : jproc=proc_dist(jterm)
553 : !Swap terms:
554 0 : call gaussfit_mpi_add_item(jterm,proc_load(iproc))
555 0 : call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
556 0 : call gaussfit_mpi_add_item(iterm,proc_load(jproc))
557 0 : call gaussfit_mpi_remove_item(jterm,proc_load(jproc))
558 : !Calculate final state
559 0 : call gaussfit_mpi_calc_deviation(deviation2,nproc,proc_load)
560 : !Swap them only if final state is better than the initial one
561 0 : if(deviation2<deviation1) then
562 0 : proc_dist(iterm)=jproc
563 0 : proc_dist(jterm)=iproc
564 : else
565 : ! Return work load to initial state
566 0 : call gaussfit_mpi_add_item(iterm,proc_load(iproc))
567 0 : call gaussfit_mpi_remove_item(jterm,proc_load(iproc))
568 0 : call gaussfit_mpi_add_item(jterm,proc_load(jproc))
569 0 : call gaussfit_mpi_remove_item(iterm,proc_load(jproc))
570 : ! write(*,*)'Back initial state'
571 : ! write(*,*)'proc_load',proc_load(:)
572 : end if
573 :
574 0 : end subroutine gaussfit_mpi_swap
575 : !!***
576 :
577 : !----------------------------------------------------------------------
578 :
579 : !!****f* m_paw_gaussfit/gaussfit_mpi_assign
580 : !! NAME
581 : !! gaussfit_mpi_assign
582 : !!
583 : !! FUNCTION
584 : !! Set task to a processor
585 : !!
586 : !! INPUTS
587 : !!
588 : !! OUTPUT
589 : !!
590 : !! SOURCE
591 :
592 0 : subroutine gaussfit_mpi_assign(iterm,nproc,nterm_bounds,&
593 0 : & proc_dist,proc_load)
594 :
595 : !Arguments ------------------------------------
596 : integer,intent(in)::iterm,nproc,nterm_bounds(2)
597 : integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc)
598 :
599 : !Local variables ------------------------------
600 : integer:: iproc,jproc,dev
601 0 : integer:: deviation(nproc),mindev
602 : character(len=100) :: msg
603 :
604 : !************************************************************************
605 :
606 0 : do iproc=1,nproc
607 : !add this term to iproc
608 0 : call gaussfit_mpi_add_item(iterm,proc_load(iproc))
609 : !calculate the deviation for this configuration
610 0 : call gaussfit_mpi_calc_deviation(dev,nproc,proc_load)
611 0 : deviation(iproc)=dev
612 : !remove this term from iproc
613 0 : call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
614 : end do
615 :
616 : !assign to jproc, the proc with minimal deviation above:
617 : jproc=-1; mindev=999999999
618 0 : do iproc=1,nproc
619 0 : if(deviation(iproc)<mindev) then
620 0 : mindev=deviation(iproc)
621 0 : jproc=iproc
622 : end if
623 : end do
624 0 : if(jproc==-1) then
625 : ! One should not get here!
626 0 : msg = 'error in accomodate_mpi'
627 0 : LIBPAW_BUG(msg)
628 : end if
629 :
630 : !assign this term for jproc
631 0 : proc_dist(iterm)=jproc
632 0 : call gaussfit_mpi_add_item(iterm,proc_load(jproc))
633 :
634 0 : end subroutine gaussfit_mpi_assign
635 : !!***
636 :
637 : !----------------------------------------------------------------------
638 :
639 : !!****f* m_paw_gaussfit/gaussfit_mpi_main
640 : !! NAME
641 : !! gaussfit_mpi_main
642 : !!
643 : !! FUNCTION
644 : !! Set charge for each processor
645 : !!
646 : !! INPUTS
647 : !!
648 : !! OUTPUT
649 : !!
650 : !! SOURCE
651 :
652 0 : subroutine gaussfit_mpi_main(nproc,nterm_bounds,proc_dist)
653 :
654 : !Arguments ------------------------------------
655 : integer,intent(in)::nproc,nterm_bounds(2)
656 : integer,intent(out)::proc_dist(nterm_bounds(1):nterm_bounds(2))
657 :
658 : !Local variables ------------------------------
659 : integer:: dev1,dev2,ii
660 : integer:: iproc,iterm,jterm,ngauss,weight
661 0 : integer:: proc_load(nproc)
662 : character(len=500) :: msg
663 :
664 : !************************************************************************
665 :
666 0 : proc_load=0; proc_dist=0 !initializations
667 :
668 : !1) get a first-trial distribution:
669 0 : do iterm=nterm_bounds(2),nterm_bounds(1),-1
670 0 : call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load)
671 : end do
672 :
673 : !Do the following 20 times
674 0 : do ii=1,20
675 : ! Calculate initial state
676 0 : call gaussfit_mpi_calc_deviation(dev1,nproc,proc_load)
677 : ! Try to swap tasks between two processors:
678 0 : do iterm=nterm_bounds(2),nterm_bounds(1),-1
679 0 : do jterm=nterm_bounds(2),nterm_bounds(1),-1
680 : call gaussfit_mpi_swap(iterm,jterm,&
681 0 : & nproc,nterm_bounds,proc_dist,proc_load)
682 : end do
683 : end do
684 : !! Try to reassign tasks to different processors
685 0 : do iterm=nterm_bounds(2),nterm_bounds(1),-1
686 0 : iproc=proc_dist(iterm)
687 : ! Remove this job from this node
688 0 : call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
689 : ! Accomodate this again:
690 0 : call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load)
691 : end do
692 : ! Calculate final state
693 0 : call gaussfit_mpi_calc_deviation(dev2,nproc,proc_load)
694 : ! If final state equals initial state, exit:
695 0 : if(dev2 == dev1) exit
696 : ! write(*,'(a)')'Deviation: ',dev2
697 : end do
698 :
699 : !Write down distribution:
700 0 : write(msg,'(3a)') 'MPI distribution',ch10,'N. gauss, iproc, weight '
701 0 : call wrtout(std_out,msg,'COLL')
702 0 : do iterm=nterm_bounds(2),nterm_bounds(1),-1
703 0 : ngauss=iterm*2
704 0 : call gaussfit_mpi_set_weight(weight,iterm)
705 0 : write(msg,'(3(i4,1x))') ngauss,proc_dist(iterm),weight
706 0 : call wrtout(std_out,msg,'COLL')
707 : end do
708 0 : write(msg,'(a)') 'Load per processor: '
709 0 : call wrtout(std_out,msg,'COLL')
710 0 : do iproc=1,nproc
711 0 : write(msg,'(i5,1x,i10)') iproc,proc_load(iproc)
712 0 : call wrtout(std_out,msg,'COLL')
713 : end do
714 :
715 0 : end subroutine gaussfit_mpi_main
716 : !!***
717 :
718 : !----------------------------------------------------------------------
719 :
720 : !!****f* m_paw_gaussfit/gaussfit_fit
721 : !! NAME
722 : !! gaussfit_fit
723 : !!
724 : !! FUNCTION
725 : !! Fits a given input function f(r) to a sum of gaussians
726 : !!
727 : !! INPUTS
728 : !! chisq= It measures how good is the fitting.
729 : !! It is defined here as sum_i^N_i f(x_i)-y(x_i)/N_i
730 : !! constrains= constraints for Gaussians
731 : !! limit(nparam)= it limits the widths of Gaussians
732 : !! maxiter=maximum number of iterations
733 : !! nparam= number of parameters (a constant times nterm)
734 : !! if(option==1)nparam=nterm*4
735 : !! if(option==2)nparam=nterm*6
736 : !! if(option==3)nparam=nterm*2
737 : !! if(option==4)nparam=nterm*4
738 : !! nterm= number of Gaussians
739 : !! nx=number of points along the x axis
740 : !! option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
741 : !! 2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
742 : !! 3 fit to a1 cos (k x^2) + a2 sin (k x^2)
743 : !! 4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
744 : !! outfile= filename for output (only written if verbosity>1)
745 : !! verbosity= controls output volume
746 : !! weight(nparam)= weights for the fitting procedure
747 : !! x(nx)= points along the x axis
748 : !! y(nx)= function to be fitted
749 : !! rpaw ,optional= paw radius
750 : !!
751 : !! OUTPUT
752 : !! y_out(nx)= fitted function
753 : !!
754 : !! SIDE EFFECTS
755 : !! if(verbosity>1) output files are written with y(x) and y_out(x)
756 : !!
757 : !! SOURCE
758 :
759 0 : subroutine gaussfit_fit(chisq,constrains,&
760 0 : & limit,maxiter,nparam,nterm,nx,option,outfile,param,&
761 0 : & verbosity,weight,x,y,y_out)
762 :
763 : !Arguments ------------------------------------
764 : integer, intent(in) :: maxiter,nparam
765 : integer,intent(in) :: nterm,nx,option,verbosity
766 : !real(dp),optional,intent(in)::rpaw
767 : !arrays
768 : integer,intent(in)::constrains(nparam)
769 : real(dp),intent(in)::limit(nparam),weight(nparam)
770 : real(dp),intent(in)::x(nx),y(nx)
771 : real(dp),intent(inout)::param(nparam)
772 : real(dp),intent(out)::chisq,y_out(nx)
773 : character(80),intent(in)::outfile
774 :
775 : !Local variables ------------------------------
776 : integer, parameter :: wfn_unit=1007
777 : integer::ix
778 : real(dp)::rerror
779 0 : real(dp),allocatable::sy(:)
780 :
781 : ! *************************************************************************
782 :
783 0 : LIBPAW_ALLOCATE(sy,(nx))
784 :
785 0 : sy(:)=1.0d0
786 :
787 : call gaussfit_rlsf(&
788 : & chisq,constrains,limit,maxiter,&
789 : & nterm,nparam,nx,option,param(1:nparam),&
790 0 : & verbosity,weight,x,y)
791 : !
792 0 : if(verbosity>=1) then
793 0 : if(option==1) then
794 0 : call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,param,x,y_out)
795 0 : elseif(option==2) then
796 0 : call gaussfit_calc_deriv_c(nparam,nterm,nx,1,param,x,y_out)
797 0 : elseif(option==3) then
798 0 : call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,param,x,y_out)
799 0 : elseif(option==4) then
800 0 : call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,param,x,y_out)
801 : end if
802 : !
803 : !
804 0 : open(wfn_unit,file=outfile,form='formatted',status='unknown')
805 : ! per_error=0.d0
806 0 : do ix=1, nx
807 0 : rerror=abs(y(ix)-y_out(ix))
808 0 : write(wfn_unit,'(6(e20.12,1x))')x(ix),y(ix),y_out(ix),rerror
809 : end do
810 0 : close(wfn_unit)
811 :
812 : end if
813 :
814 0 : LIBPAW_DEALLOCATE(sy)
815 :
816 0 : end subroutine gaussfit_fit
817 : !!***
818 :
819 : !----------------------------------------------------------------------
820 :
821 : !!****f* m_paw_gaussfit/gaussfit_calc_deriv_r
822 : !! NAME
823 : !! gaussfit_calc_deriv_r
824 : !!
825 : !! FUNCTION
826 : !! Calculate derivatives for Gaussians
827 : !! Only relevant for fitting Gaussians algorithm.
828 : !! The Gaussians expressions are defined in the comments of "gaussfit_main"
829 : !!
830 : !! INPUTS
831 : !!
832 : !! OUTPUT
833 : !!
834 : !! SOURCE
835 :
836 : subroutine gaussfit_calc_deriv_r(nterm,nparam,nx,opt,param,x,y_out,&
837 : & deriv) ! optional
838 :
839 : !Arguments -------------------------------
840 : integer,intent(in)::nx !number of point in the x grid
841 : integer,intent(in)::nparam !number of parameters
842 : integer,intent(in)::nterm !number of gaussian expressions
843 : integer,intent(in)::opt !option:
844 : !1) calculate only f(x)
845 : !2) calculate f(x) and its derivatives
846 : real(dp),intent(in)::param(nparam) !parameters
847 : real(dp),intent(in)::x(nx) !xgrid
848 : real(dp),intent(out)::y_out(nx) !f(x)
849 : real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
850 :
851 : !Local variables-------------------------------
852 : integer::iexp,ii
853 : real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
854 : real(dp)::term1(nx,nterm)
855 : real(dp)::aux1(nx)
856 : !real(dp)::step
857 :
858 : ! *********************************************************************
859 :
860 : !
861 : !Initialize
862 : !
863 : y_out(:)=0.d0
864 : !
865 : !Get parameters from parameters array:
866 : !
867 : alpha1(:)=param(1:nterm)
868 : alpha2(:)=param(nterm+1:2*nterm)
869 : alpha3(:)=param(2*nterm+1:3*nterm)
870 : !
871 : !alpha3
872 : !set to constant values of x
873 : !step=rpaw/real(nterm,dp)
874 : !do ii=1,nterm
875 : !raux=step*real(ii,dp)
876 : !alpha3(ii)=raux
877 : !end do
878 : !
879 : !
880 : !
881 : !calculate useful quantities
882 : !
883 : do iexp=1,nterm
884 : aux1(:)=-alpha2(iexp)*(x(:)-alpha3(iexp))**2
885 : term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
886 : end do
887 : !
888 : do iexp=1,nterm
889 : y_out(:)=y_out(:)+term1(:,iexp)
890 : end do
891 : !
892 : !Calculate derivatives:
893 : !
894 : if(opt==2) then
895 : !
896 : ! alpha1
897 : !
898 : do iexp=1,nterm
899 : aux1(:)=term1(:,iexp)/alpha1(iexp)
900 : deriv(:,iexp)=aux1(:)
901 : end do
902 : !
903 : ! alpha2
904 : !
905 : do iexp=1,nterm
906 : ii=nterm+iexp
907 : aux1(:)=-term1(:,iexp)*(x(:)-alpha3(iexp))
908 : deriv(:,ii)=aux1(:)
909 : deriv(:,ii)=0.1d0
910 : end do
911 : !
912 : ! alpha3
913 : !
914 : do iexp=1,nterm
915 : ii=2*nterm+iexp
916 : aux1(:)=term1(:,iexp)*2.d0*alpha2(iexp)
917 : aux1(:)=aux1(:)*(x(:)-alpha3(iexp))
918 : deriv(:,ii)=aux1(:)
919 : end do
920 : end if
921 :
922 : end subroutine gaussfit_calc_deriv_r
923 : !!***
924 :
925 : !----------------------------------------------------------------------
926 :
927 : !!****f* m_paw_gaussfit/gaussfit_calc_deriv_c3
928 : !! NAME
929 : !! gaussfit_calc_deriv_c3
930 : !!
931 : !! FUNCTION
932 : !! Calculate expressions and derivatives for Gaussians fitting.
933 : !! The Gaussians expressions are defined in the comments of "gaussfit_main"
934 : !!
935 : !! INPUTS
936 : !!
937 : !! OUTPUT
938 : !!
939 : !! SOURCE
940 :
941 0 : subroutine gaussfit_calc_deriv_c3(nparam,nterm,nx,opt,param,x,y_out,&
942 0 : & deriv) ! optional
943 :
944 : !Arguments -------------------------------
945 : integer,intent(in)::nparam !number of parameters
946 : integer,intent(in)::nterm !number of gaussian expressions
947 : integer,intent(in)::nx !number of point in the x grid
948 : integer,intent(in)::opt !option:
949 : !1) calculate only f(x)
950 : !2) calculate f(x) and its derivatives
951 : real(dp),intent(in)::param(nparam) !parameters
952 : real(dp),intent(in)::x(nx) !xgrid
953 : real(dp),intent(out)::y_out(nx) !f(x)
954 : real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
955 :
956 : !Local variables-------------------------------
957 : integer::iexp,ii
958 : real(dp)::sep
959 0 : real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
960 0 : real(dp)::term1(nx,nterm)
961 0 : real(dp)::sin1(nx,nterm),cos1(nx,nterm)
962 0 : real(dp)::aux1(nx),aux2(nx)
963 :
964 : ! *********************************************************************
965 :
966 : !
967 : !Initialize
968 : !
969 0 : y_out(:)=0.d0
970 : !
971 : sep=1.2d0
972 : !
973 : !Get param from param array:
974 : !
975 0 : alpha1(:)=param(1:nterm)
976 0 : alpha2(:)=param(nterm+1:2*nterm)
977 : !
978 0 : do ii=1,nterm
979 0 : alpha3(ii)=sep**(ii)
980 : end do
981 : !
982 : !calculate useful quantities
983 : !
984 0 : do iexp=1,nterm
985 0 : aux1(:)=alpha3(iexp)*x(:)**2
986 : !
987 0 : sin1(:,iexp)=sin(aux1(:))
988 0 : cos1(:,iexp)=cos(aux1(:))
989 : end do
990 : !
991 0 : do iexp=1,nterm
992 0 : aux1(:)=alpha1(iexp)*sin1(:,iexp)
993 0 : aux2(:)=alpha2(iexp)*cos1(:,iexp)
994 0 : term1(:,iexp)=aux1(:)+aux2(:)
995 0 : y_out(:)=y_out(:)+term1(:,iexp)
996 : end do
997 : !
998 : !Calculate derivatives:
999 : !
1000 0 : if(opt==2) then
1001 : !
1002 : ! alpha1
1003 : !
1004 0 : do iexp=1,nterm
1005 0 : deriv(:,iexp)=sin1(:,iexp)
1006 : end do
1007 : !
1008 : ! alpha2
1009 : !
1010 0 : do iexp=1,nterm
1011 0 : ii=nterm+iexp
1012 0 : deriv(:,ii)=cos1(:,iexp)
1013 : end do
1014 : end if
1015 :
1016 0 : end subroutine gaussfit_calc_deriv_c3
1017 : !!***
1018 :
1019 : !----------------------------------------------------------------------
1020 :
1021 : !!****f* m_paw_gaussfit/gaussfit_calc_deriv_c2
1022 : !! NAME
1023 : !! gaussfit_calc_deriv_c2
1024 : !!
1025 : !! FUNCTION
1026 : !! Calculate expressions and derivatives for Gaussians fitting.
1027 : !! The Gaussians expressions are defined in the comments of "gaussfit_main"
1028 : !!
1029 : !! INPUTS
1030 : !!
1031 : !! OUTPUT
1032 : !!
1033 : !! SOURCE
1034 :
1035 0 : subroutine gaussfit_calc_deriv_c2(nparam,nterm,nx,opt,param,x,y_out,&
1036 0 : & deriv) ! optional
1037 :
1038 : !Arguments -------------------------------
1039 : integer,intent(in)::nparam !number of param
1040 : integer,intent(in)::nterm !number of gaussian expressions
1041 : integer,intent(in)::nx !number of point in the x grid
1042 : integer,intent(in)::opt !option:
1043 : !1) calculate only f(x)
1044 : !2) calculate f(x) and its derivatives
1045 : real(dp),intent(in)::param(nparam) !parameters
1046 : real(dp),intent(in)::x(nx) !xgrid
1047 : real(dp),intent(out)::y_out(nx) !f(x)
1048 : real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1049 :
1050 : !Local variables-------------------------------
1051 : integer::iexp,ii
1052 0 : real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1053 0 : real(dp)::alpha4(nterm)
1054 0 : real(dp)::term1(nx,nterm)
1055 0 : real(dp)::sin1(nx,nterm),sin2(nx,nterm),cos1(nx,nterm),cos2(nx,nterm)
1056 0 : real(dp)::aux1(nx),aux2(nx)
1057 :
1058 : ! *********************************************************************
1059 :
1060 : !
1061 : !Initialize
1062 : !
1063 0 : y_out(:)=0.d0
1064 : !
1065 : !Get param from param array:
1066 : !
1067 0 : alpha1(:)=param(1:nterm)
1068 0 : alpha2(:)=param(nterm+1:2*nterm)
1069 0 : alpha3(:)=param(2*nterm+1:3*nterm)
1070 0 : alpha4(:)=param(3*nterm+1:4*nterm)
1071 : !
1072 : !calculate useful quantities
1073 : !
1074 : !
1075 0 : do iexp=1,nterm
1076 0 : aux1(:)=alpha2(iexp)*x(:)**2
1077 0 : sin1(:,iexp)=sin(aux1(:))
1078 : !
1079 0 : aux1(:)=alpha4(iexp)*x(:)**2
1080 0 : sin2(:,iexp)=sin(aux1(:))
1081 : !
1082 0 : aux1(:)=alpha2(iexp)*x(:)**2
1083 0 : cos1(:,iexp)=cos(aux1(:))
1084 : !
1085 0 : aux1(:)=alpha4(iexp)*x(:)**2
1086 0 : cos2(:,iexp)=cos(aux1(:))
1087 : end do
1088 : !
1089 0 : do iexp=1,nterm
1090 0 : aux1(:)=alpha1(iexp)*sin1(:,iexp)
1091 0 : aux2(:)=alpha3(iexp)*cos2(:,iexp)
1092 0 : term1(:,iexp)=aux1(:)+aux2(:)
1093 0 : y_out(:)=y_out(:)+term1(:,iexp)
1094 : end do
1095 : !
1096 : !Calculate derivatives:
1097 : !
1098 0 : if(opt==2) then
1099 : !
1100 : ! alpha1
1101 : !
1102 0 : do iexp=1,nterm
1103 0 : deriv(:,iexp)=sin1(:,iexp)
1104 : end do
1105 : !
1106 : ! alpha2
1107 : !
1108 0 : do iexp=1,nterm
1109 0 : ii=nterm+iexp
1110 0 : aux1(:)=alpha1(iexp)*cos1(:,iexp)*x(:)**2
1111 0 : deriv(:,ii)=aux1(:)
1112 : end do
1113 : !
1114 : ! alpha3
1115 : !
1116 0 : do iexp=1,nterm
1117 0 : ii=2*nterm+iexp
1118 0 : deriv(:,ii)=cos2(:,iexp)
1119 : end do
1120 : !
1121 : ! alpha4
1122 : !
1123 0 : do iexp=1,nterm
1124 0 : ii=3*nterm+iexp
1125 0 : aux1(:)=-alpha3(iexp)*sin2(:,iexp)*x(:)**2
1126 0 : deriv(:,ii)=aux1(:)
1127 : end do
1128 : end if
1129 :
1130 0 : end subroutine gaussfit_calc_deriv_c2
1131 : !!***
1132 :
1133 : !----------------------------------------------------------------------
1134 :
1135 : !!****f* m_paw_gaussfit/gaussfit_calc_deriv_c
1136 : !! NAME
1137 : !! gaussfit_calc_deriv_c
1138 : !!
1139 : !! FUNCTION
1140 : !! Calculate expressions and derivatives for Gaussians fitting.
1141 : !! The Gaussians expressions are defined in the comments of "gaussfit_main"
1142 : !!
1143 : !! INPUTS
1144 : !!
1145 : !! OUTPUT
1146 : !!
1147 : !! SOURCE
1148 :
1149 0 : subroutine gaussfit_calc_deriv_c(nparam,nterm,nx,opt,param,x,y_out,&
1150 0 : & deriv) ! optional
1151 :
1152 : !Arguments -------------------------------
1153 : integer,intent(in)::nparam !number of parameters
1154 : integer,intent(in)::nterm !number of gaussian expressions
1155 : integer,intent(in)::nx !number of point in the x grid
1156 : integer,intent(in)::opt !option:
1157 : !1) calculate only f(x)
1158 : !2) calculate f(x) and its derivatives
1159 : real(dp),intent(in)::param(nparam) !parameters
1160 : real(dp),intent(in)::x(nx) !xgrid
1161 : real(dp),intent(out)::y_out(nx) !f(x)
1162 : real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1163 :
1164 : !Local variables-------------------------------
1165 : integer::iexp,ii
1166 0 : real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1167 0 : real(dp)::alpha4(nterm),alpha5(nterm),alpha6(nterm)
1168 0 : real(dp)::aux1(nx),aux2(nx)
1169 0 : real(dp)::cos1(nx,nterm),cos2(nx,nterm),sin1(nx,nterm),sin2(nx,nterm)
1170 0 : real(dp)::term1(nx,nterm),term2(nx,nterm)
1171 :
1172 : ! *********************************************************************
1173 :
1174 : !
1175 : !Initialize
1176 : !
1177 0 : y_out(:)=0.d0
1178 : !
1179 : !Get parameters from param array:
1180 : !
1181 0 : alpha1(:)=param(1:nterm)
1182 0 : alpha2(:)=param(nterm+1:2*nterm)
1183 0 : alpha3(:)=param(2*nterm+1:3*nterm)
1184 0 : alpha4(:)=param(3*nterm+1:4*nterm)
1185 0 : alpha5(:)=param(4*nterm+1:5*nterm)
1186 0 : alpha6(:)=param(5*nterm+1:6*nterm)
1187 : !
1188 : !calculate useful quantities
1189 : !
1190 0 : do iexp=1,nterm
1191 0 : aux1(:)=-alpha2(iexp)*x(:)**2
1192 0 : term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
1193 : end do
1194 : !
1195 0 : do iexp=1,nterm
1196 0 : aux1(:)=alpha4(iexp)*x(:)**2
1197 0 : sin1(:,iexp)=sin(aux1(:))
1198 : !
1199 0 : aux1(:)=alpha6(iexp)*x(:)**2
1200 0 : sin2(:,iexp)=sin(aux1(:))
1201 : !
1202 0 : aux1(:)=alpha4(iexp)*x(:)**2
1203 0 : cos1(:,iexp)=cos(aux1(:))
1204 : !
1205 0 : aux1(:)=alpha6(iexp)*x(:)**2
1206 0 : cos2(:,iexp)=cos(aux1(:))
1207 : end do
1208 : !
1209 0 : do iexp=1,nterm
1210 0 : aux1(:)=alpha3(iexp)*sin1(:,iexp)
1211 0 : aux2(:)=alpha5(iexp)*cos2(:,iexp)
1212 0 : term2(:,iexp)=aux1(:)+aux2(:)
1213 0 : y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp)
1214 : end do
1215 : !
1216 : !Calculate derivatives:
1217 : !
1218 0 : if(opt==2) then
1219 : !
1220 : ! alpha1
1221 : !
1222 0 : do iexp=1,nterm
1223 0 : aux1(:)=term1(:,iexp)/alpha1(iexp)
1224 0 : aux2(:)=aux1(:)*term2(:,iexp)
1225 0 : deriv(:,iexp)=aux2(:)
1226 : end do
1227 : !
1228 : ! alpha2
1229 : !
1230 0 : do iexp=1,nterm
1231 0 : ii=nterm+iexp
1232 0 : aux1(:)=-term1(:,iexp)*term2(:,iexp)
1233 0 : aux2(:)=aux1(:)*x(:)**2
1234 0 : deriv(:,ii)=aux2(:)
1235 : end do
1236 : !
1237 : ! alpha3
1238 : !
1239 0 : do iexp=1,nterm
1240 0 : ii=2*nterm+iexp
1241 0 : aux1(:)=term1(:,iexp)*sin1(:,iexp)
1242 0 : deriv(:,ii)=aux1(:)
1243 : end do
1244 : !
1245 : ! alpha4
1246 : !
1247 0 : do iexp=1,nterm
1248 0 : ii=3*nterm+iexp
1249 0 : aux1(:)=term1(:,iexp)*alpha3(iexp)
1250 0 : aux2(:)=cos1(:,iexp)*x(:)**2
1251 0 : deriv(:,ii)=aux2(:)*aux1(:)
1252 : end do
1253 : !
1254 : ! alpha5
1255 : !
1256 0 : do iexp=1,nterm
1257 0 : ii=4*nterm+iexp
1258 0 : aux1(:)=term1(:,iexp)*cos2(:,iexp)
1259 0 : deriv(:,ii)=aux1(:)
1260 : end do
1261 : !
1262 : ! alpha6
1263 : !
1264 0 : do iexp=1,nterm
1265 0 : ii=5*nterm+iexp
1266 0 : aux1(:)=-term1(:,iexp)*alpha5(iexp)
1267 0 : aux2(:)=sin2(:,iexp)*x(:)**2
1268 0 : deriv(:,ii)=aux1(:)*aux2(:)
1269 : end do
1270 : end if
1271 :
1272 0 : end subroutine gaussfit_calc_deriv_c
1273 : !!***
1274 :
1275 : !----------------------------------------------------------------------
1276 :
1277 : !!****f* m_paw_gaussfit/gaussfit_calc_deriv_c4
1278 : !! NAME
1279 : !! gaussfit_calc_deriv_c4
1280 : !!
1281 : !! FUNCTION
1282 : !! Calculate expressions and derivatives for Gaussians fitting.
1283 : !! The Gaussians expressions are defined in the comments of "gaussfit_main"
1284 : !!
1285 : !! INPUTS
1286 : !!
1287 : !! OUTPUT
1288 : !!
1289 : !! SOURCE
1290 :
1291 0 : subroutine gaussfit_calc_deriv_c4(nparam,nterm,nx,opt,param,x,y_out,&
1292 0 : & deriv) ! optional
1293 :
1294 : !Arguments -------------------------------
1295 : integer,intent(in)::nparam !number of parameters
1296 : integer,intent(in)::nterm !number of gaussian expressions
1297 : integer,intent(in)::nx !number of point in the x grid
1298 : integer,intent(in)::opt !option:
1299 : !1) calculate only f(x)
1300 : !2) calculate f(x) and its derivatives
1301 : real(dp),intent(in)::param(nparam) !parameters
1302 : real(dp),intent(in)::x(nx) !xgrid
1303 : real(dp),intent(out)::y_out(nx) !f(x)
1304 : real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1305 :
1306 : !Local variables-------------------------------
1307 : integer::iexp,ii
1308 : real(dp)::raux,sep
1309 0 : real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1310 0 : real(dp)::alpha4(nterm),alpha5(nterm)
1311 0 : real(dp)::aux1(nx),aux2(nx)
1312 0 : real(dp)::cos1(nx,nterm),sin1(nx,nterm)
1313 0 : real(dp)::term1(nx,nterm),term2(nx,nterm)
1314 :
1315 : ! *********************************************************************
1316 :
1317 : !
1318 : !Initialize
1319 : !
1320 0 : sep=1.1d0
1321 0 : y_out(:)=0.d0
1322 :
1323 : !Get parameters from param array:
1324 : !
1325 0 : alpha1(:)=param(1:nterm)
1326 0 : alpha2(:)=param(nterm+1:2*nterm)
1327 0 : alpha3(:)=param(2*nterm+1:3*nterm)
1328 0 : alpha4(:)=param(3*nterm+1:4*nterm)
1329 : !
1330 : !
1331 0 : raux=(2.d0*pi)/real(nterm,dp)
1332 0 : do ii=1,nterm
1333 0 : alpha5(ii)=sep**(ii)
1334 : ! alpha5(ii)=raux*real(ii-1,dp)
1335 : end do
1336 : !
1337 : !calculate useful quantities
1338 : !
1339 0 : do iexp=1,nterm
1340 0 : aux1(:)=-alpha2(iexp)*x(:)**2
1341 0 : term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
1342 : end do
1343 : !
1344 0 : do iexp=1,nterm
1345 0 : aux1(:)=alpha5(iexp)*x(:)**2
1346 : !
1347 0 : sin1(:,iexp)=sin(aux1(:))
1348 0 : cos1(:,iexp)=cos(aux1(:))
1349 : end do
1350 : !
1351 0 : do iexp=1,nterm
1352 0 : aux1(:)=alpha3(iexp)*sin1(:,iexp)
1353 0 : aux2(:)=alpha4(iexp)*cos1(:,iexp)
1354 0 : term2(:,iexp)=aux1(:)+aux2(:)
1355 0 : y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp)
1356 : end do
1357 : !
1358 : !Calculate derivatives:
1359 : !
1360 0 : if(opt==2) then
1361 : !
1362 : ! alpha1
1363 : !
1364 0 : do iexp=1,nterm
1365 0 : aux1(:)=term1(:,iexp)/alpha1(iexp)
1366 0 : aux2(:)=aux1(:)*term2(:,iexp)
1367 0 : deriv(:,iexp)=aux2(:)
1368 : end do
1369 : !
1370 : ! alpha2
1371 : !
1372 0 : do iexp=1,nterm
1373 0 : ii=nterm+iexp
1374 0 : aux1(:)=-term1(:,iexp)*term2(:,iexp)
1375 0 : aux2(:)=aux1(:)*x(:)**2
1376 0 : deriv(:,ii)=aux2(:)
1377 : end do
1378 : !
1379 : ! alpha3
1380 : !
1381 0 : do iexp=1,nterm
1382 0 : ii=2*nterm+iexp
1383 0 : aux1(:)=term1(:,iexp)*sin1(:,iexp)
1384 0 : deriv(:,ii)=aux1(:)
1385 : end do
1386 : !
1387 : ! alpha4
1388 : !
1389 0 : do iexp=1,nterm
1390 0 : ii=3*nterm+iexp
1391 0 : aux1(:)=term1(:,iexp)*cos1(:,iexp)
1392 0 : deriv(:,ii)=aux1(:)
1393 : end do
1394 : end if
1395 :
1396 0 : end subroutine gaussfit_calc_deriv_c4
1397 : !!***
1398 :
1399 : !----------------------------------------------------------------------
1400 :
1401 : !!****f* m_paw_gaussfit/gaussfit_rlsf
1402 : !! NAME
1403 : !! gaussfit_rlsf
1404 : !!
1405 : !! FUNCTION
1406 : !! Fits a given function to a sum of Gaussians.
1407 : !! Uses the Levenberg-Marquardt algorithm.
1408 : !!
1409 : !! COPYRIGHT
1410 : !! Copyright (C) 2011-2026 ABINIT group (T. Rangel)
1411 : !! This file is distributed under the terms of the
1412 : !! GNU General Public License, see ~abinit/COPYING
1413 : !! or http://www.gnu.org/copyleft/gpl.txt .
1414 : !! The original Levemberg Marquardt routines were written by Armando Sole
1415 : !! These were modified for the ARPUS spectra in the BigDFT code by A. Mirone.
1416 : !! These were re-writen in Fortran and further modified in ABINIT for our particular needs.
1417 : !!
1418 : !! INPUTS
1419 : !! option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
1420 : !! 2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
1421 : !! 3 fit to a1 cos (k x^2) + a2 sin (k x^2)
1422 : !! 4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
1423 : !! if(option==1)mparam=nterm_bounds(2)*4
1424 : !! if(option==2)mparam=nterm_bounds(2)*6
1425 : !! if(option==3)mparam=nterm_bounds(2)*2
1426 : !! if(option==4)mparam=nterm_bounds(2)*4
1427 : !!
1428 : !! OUTPUT
1429 : !!
1430 : !! SOURCE
1431 :
1432 0 : subroutine gaussfit_rlsf(&
1433 0 : &chisq,constrains,limit,maxiter,&
1434 0 : &nterm,nparam,nx,option,parameters,&
1435 0 : &verbosity,weight,x,y)
1436 :
1437 : !Arguments -------------------------------
1438 : real(dp),parameter::deltachi=tol10
1439 : integer, intent(in) ::maxiter,nparam,nterm,nx
1440 : integer, intent(in) ::option ,verbosity
1441 : integer, intent(in) ::constrains(nparam)
1442 : real(dp),intent(out)::chisq
1443 : !arrays
1444 : real(dp),intent(in)::limit(nparam),weight(nparam)
1445 : real(dp),intent(inout)::parameters(nparam)
1446 : real(dp),intent(in)::x(nx),y(nx)
1447 :
1448 : !Local variables-------------------------------
1449 : integer::flag,ii,info,iter,jj,niter
1450 : real(dp):: deltax
1451 : real(dp)::chisq0,flambda,eta,lastdeltachi
1452 0 : integer::ipvt(nparam)
1453 0 : real(dp)::alpha(nparam,nparam)
1454 0 : real(dp)::alpha0(nparam,nparam),beta(nparam)
1455 0 : real(dp)::deltapar(nparam)
1456 0 : real(dp)::tmp1(nparam,nparam)
1457 0 : real(dp)::work(nparam)
1458 0 : real(dp)::workpar(nparam)
1459 0 : real(dp)::yfit(nx)
1460 : character(len=500) :: msg
1461 :
1462 : ! *********************************************************************
1463 :
1464 : !
1465 : !flambda=1e-6
1466 0 : flambda=1e-7
1467 : !iter=maxiter !later it is changed
1468 0 : niter=0
1469 0 : deltax=x(2)-x(1) !we assume this is a linear grid
1470 : !
1471 0 : iter_loop: do iter=1,maxiter
1472 : !
1473 : call gaussfit_chisq_alpha_beta(alpha0,beta,chisq0,&
1474 0 : & nparam,nterm,nx,option,parameters,x,y)
1475 : !
1476 0 : flag=0
1477 0 : lastdeltachi=chisq0
1478 : !
1479 : !
1480 0 : while_flag: do
1481 : if(flag .ne. 0) exit while_flag
1482 : !
1483 0 : tmp1=0.d0
1484 0 : do ii=1,nparam
1485 0 : tmp1(ii,ii)=1.d0*flambda !identity matrix * flambda
1486 : end do
1487 0 : alpha=alpha0+tmp1*alpha0
1488 : ! Invert alpha matrix
1489 0 : tmp1=alpha
1490 0 : call dgetrf(nparam,nparam,tmp1,nparam,ipvt,info)
1491 0 : if (.not.info==0) then
1492 0 : if(verbosity>1) then
1493 0 : write(msg,'(a)')'Matrix is singular'
1494 0 : call wrtout(std_out,msg,'COLL')
1495 : end if
1496 0 : chisq=-1.d0
1497 0 : exit iter_loop
1498 : end if
1499 0 : call dgetri(nparam,tmp1,nparam,ipvt,work,nparam,info)
1500 0 : deltapar=0.d0
1501 0 : if (.not.info==0) then
1502 0 : if(verbosity>2) then
1503 0 : write(msg,'(a)')'Matrix is singular'
1504 0 : call wrtout(std_out,msg,'COLL')
1505 : end if
1506 0 : chisq=-1.d0
1507 0 : exit iter_loop
1508 : end if
1509 : !
1510 0 : if(tmp1(1,1) .ne. tmp1(1,1)) then !If is NaN
1511 0 : chisq=-1.d0
1512 0 : exit iter_loop
1513 : end if
1514 0 : if(abs(tmp1(1,1)) == tmp1(1,1)*tmp1(1,1)) then !If is infinity
1515 0 : chisq=-1.d0
1516 0 : exit iter_loop
1517 : end if
1518 : !
1519 0 : do ii=1,nparam
1520 0 : do jj=1,nparam
1521 0 : deltapar(ii)=deltapar(ii)+beta(jj)*tmp1(jj,ii)
1522 : end do
1523 : end do
1524 : ! apply constrains
1525 0 : workpar(1:nparam)=parameters(1:nparam)+deltapar(1:nparam)*weight(1:nparam)
1526 0 : call gaussfit_apply_constrains(constrains,limit,nparam,workpar)
1527 : !
1528 0 : if(option==1) then
1529 0 : call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,workpar,x,yfit)
1530 0 : elseif(option==2) then
1531 0 : call gaussfit_calc_deriv_c(nparam,nterm,nx,1,workpar,x,yfit)
1532 0 : elseif(option==3) then
1533 0 : call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,workpar,x,yfit)
1534 0 : elseif(option==4) then
1535 0 : call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,workpar,x,yfit)
1536 : end if
1537 0 : chisq=0.d0
1538 0 : do ii=1,nx
1539 0 : chisq=chisq + ((y(ii)-yfit(ii)))**2
1540 : end do
1541 0 : chisq=chisq*deltax
1542 : !
1543 : ! write(*,'("chisq ",f12.5," chisq0 ",f12.5)')chisq,chisq0
1544 : !
1545 0 : if(chisq > chisq0) then
1546 0 : flambda=flambda*2.0d0
1547 0 : if( flambda > 1000.d0) then
1548 0 : flag=1
1549 : ! iter=0
1550 0 : if(verbosity>2) then
1551 0 : write(msg,'(a)')'flambda > 1000.d0'
1552 0 : call wrtout(std_out,msg,'COLL')
1553 : end if
1554 : exit iter_loop
1555 : end if
1556 : else
1557 0 : flag=1
1558 0 : parameters=workpar
1559 0 : eta=0.d0
1560 0 : lastdeltachi=(chisq0-chisq) !/(chisq0+eta)
1561 0 : if(lastdeltachi<deltachi) cycle
1562 0 : chisq0=chisq
1563 0 : flambda=flambda/2.d0
1564 0 : if(verbosity>2) then
1565 0 : write(msg,'("iter = ",i4," chisq = ",e15.6)')iter,chisq
1566 0 : call wrtout(std_out,msg,'COLL')
1567 : end if
1568 : end if
1569 : end do while_flag
1570 : end do iter_loop
1571 :
1572 0 : end subroutine gaussfit_rlsf
1573 : !!***
1574 :
1575 : !----------------------------------------------------------------------
1576 :
1577 : !!****f* m_paw_gaussfit/gaussfit_chisq_alpha_beta
1578 : !! NAME
1579 : !! gaussfit_chisq_alpha_beta
1580 : !!
1581 : !! FUNCTION
1582 : !! Finds chisq, alpha and beta parameters for LSF using the Levenberg-Marquardt algorithm.
1583 : !!
1584 : !! COPYRIGHT
1585 : !! Copyright (C) 2011-2026 ABINIT group (T. Rangel)
1586 : !! This file is distributed under the terms of the
1587 : !! GNU General Public License, see ~abinit/COPYING
1588 : !! or http://www.gnu.org/copyleft/gpl.txt .
1589 : !! The original Levemberg Marquardt routines were written by Armando Sole
1590 : !! These were modified for the ARPUS spectra in the BigDFT code by A. Mirone.
1591 : !! These were re-writen in Fortran and further modified in ABINIT for our particular needs.
1592 : !!
1593 : !! INPUTS
1594 : !!
1595 : !! OUTPUT
1596 : !!
1597 : !! SOURCE
1598 :
1599 0 : subroutine gaussfit_chisq_alpha_beta(alpha,beta,chisq,&
1600 0 : & nparam,nterm,nx,option,parameters,x,y)
1601 :
1602 : !Arguments -------------------------------
1603 : integer,intent(in)::nparam,nterm,nx
1604 : integer,intent(in)::option
1605 : real(dp),intent(out)::chisq
1606 : real(dp),intent(in)::parameters(nparam),x(nx),y(nx)
1607 : real(dp),intent(out)::alpha(nparam,nparam),beta(nparam)
1608 :
1609 : !Local variables-------------------------------
1610 : integer::ii,jj,kk
1611 : real(dp)::deltax,help1
1612 : !arrays
1613 0 : real(dp)::deltay(nx),deriv(nx,nparam),derivi(nx)
1614 0 : real(dp)::yfit(nx)
1615 0 : real(dp)::help0(nx),help2(nx),help3(nparam)
1616 :
1617 : ! *********************************************************************
1618 :
1619 0 : deltax=x(2)-x(1) !we assume a linear grid
1620 : !
1621 0 : if(option==1) then
1622 0 : call gaussfit_calc_deriv_c2(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1623 0 : elseif(option==2) then
1624 0 : call gaussfit_calc_deriv_c(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1625 0 : elseif(option==3) then
1626 0 : call gaussfit_calc_deriv_c3(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1627 0 : elseif(option==4) then
1628 0 : call gaussfit_calc_deriv_c4(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1629 : end if
1630 0 : deltay=y-yfit
1631 0 : help0=deltay
1632 : !
1633 0 : do ii=1,nparam
1634 0 : derivi(:)=deriv(:,ii)
1635 : help1=0.d0
1636 0 : do jj=1,nx
1637 0 : help1=help1+help0(jj)*derivi(jj)
1638 : end do
1639 0 : beta(ii)=help1
1640 : ! help1 = innerproduct(deriv,weight*derivi)
1641 : ! below I use help3 instead for the array dimenstions
1642 0 : help3=0.d0
1643 0 : do kk=1,nparam
1644 0 : do jj=1,nx
1645 0 : help3(kk)=help3(kk)+deriv(jj,kk)*derivi(jj)
1646 : end do
1647 : end do
1648 : ! !
1649 0 : alpha(:,ii)=help3(:)
1650 : end do
1651 : !
1652 0 : help2(:)=help0(:)*deltay(:)
1653 0 : chisq=sum(help2)
1654 0 : chisq=chisq*deltax
1655 :
1656 0 : end subroutine gaussfit_chisq_alpha_beta
1657 : !!***
1658 :
1659 : !----------------------------------------------------------------------
1660 :
1661 : !!****f* m_paw_gaussfit/gaussfit_set_param1
1662 : !! NAME
1663 : !! gaussfit_set_param1
1664 : !!
1665 : !! FUNCTION
1666 : !! Sets parameters for LSF
1667 : !!
1668 : !! INPUTS
1669 : !!
1670 : !! OUTPUT
1671 : !!
1672 : !! SOURCE
1673 :
1674 : subroutine gaussfit_set_param1(nterm,nparam,nx,param,sep,x,y)
1675 :
1676 : !Arguments -------------------------------
1677 : integer,intent(in)::nterm,nparam,nx
1678 : real(dp),intent(in)::sep
1679 : real(dp),intent(in)::x(nx),y(nx)
1680 : real(dp),intent(out)::param(nparam)
1681 :
1682 : !Local variables-------------------------------
1683 : integer::ii,jj
1684 : real(dp)::raux
1685 :
1686 : ! *********************************************************************
1687 :
1688 : !
1689 : param(:)=1.0d0
1690 : !exps=1.0/(x(nx)**2)
1691 : !
1692 : !alpha1
1693 :
1694 : !raux=maxval( y(:),nx )
1695 : !maxval gives problems in some architectures:
1696 : raux=-9999999
1697 : do ii=1,nx
1698 : if(raux<y(ii)) raux=y(ii)
1699 : end do
1700 : param(1:nterm)=raux
1701 : !
1702 : !alpha2
1703 : !
1704 : !y(r_c)=e^{-\alpha1 r_c}
1705 : raux=-log(abs(y(nx))+tol10)/(x(nx)**2)
1706 : param(nterm+1:nterm+2)=raux
1707 : !
1708 : !raux=0.5d0*pi/real(nterm,dp)
1709 : do jj=1,nterm
1710 : ii=jj+3*nterm
1711 : param(ii)=sep**(jj)
1712 : ! param(ii)=raux*real(jj,dp)
1713 : end do
1714 : !
1715 : do jj=1,nterm
1716 : ii=jj+5*nterm
1717 : param(ii)=sep**(jj)
1718 : ! param(ii)=raux*real(jj,dp)
1719 : end do
1720 :
1721 : end subroutine gaussfit_set_param1
1722 : !!***
1723 :
1724 : !----------------------------------------------------------------------
1725 :
1726 : !!****f* m_paw_gaussfit/gaussfit_set_param2
1727 : !! NAME
1728 : !! gaussfit_set_param2
1729 : !!
1730 : !! FUNCTION
1731 : !! Sets parameters for LSF
1732 : !!
1733 : !! INPUTS
1734 : !!
1735 : !! OUTPUT
1736 : !!
1737 : !! SOURCE
1738 :
1739 : subroutine gaussfit_set_param2(nterm,nparam,nx,param,rpaw,x,y)
1740 :
1741 : !Arguments -------------------------------
1742 : integer,intent(in)::nterm,nparam,nx
1743 : real(dp),intent(in)::rpaw
1744 : real(dp),intent(in)::x(nx),y(nx)
1745 : real(dp),intent(out)::param(nparam)
1746 :
1747 : !Local variables-------------------------------
1748 : integer::ii,jj
1749 : real(dp)::exps,raux,sig,step
1750 :
1751 : ! *************************************************************************
1752 :
1753 : step=rpaw/real(nterm-1,dp)
1754 : !exps=1.0/(rpaw/(real(nterm,dp)/2.d0))**2
1755 : exps=1.0/(step*1.0d0)**2
1756 : !alpha2 (width of gaussians)
1757 : !Set to exps*real(nterm,dp)**2, for a good guess
1758 : param(nterm+1:2*nterm)=exps !*real(nterm,dp)**2
1759 : !
1760 : do jj=1,nterm
1761 : ! alpha3
1762 : ! set to constant values of x
1763 : ii=jj+2*nterm
1764 : raux=step*real(jj-1,dp)
1765 : param(ii)=raux
1766 : ! alpha1
1767 : ! set to the value of y at that point
1768 : call gaussfit_param2_findsign()
1769 : param(jj)=sig
1770 : end do
1771 : !
1772 : contains
1773 : !!***
1774 :
1775 : !!****f* gaussfit_set_param2/gaussfit_param2_findsign
1776 : !! NAME
1777 : !! gaussfit_param2_findsign
1778 : !!
1779 : !! FUNCTION
1780 : !! Finds the value of y at a given point
1781 : !! This was taken out of gaussfit_set_param2 to make the code more
1782 : !! readable.
1783 : !!
1784 : !! INPUTS
1785 : !!
1786 : !! OUTPUT
1787 : !!
1788 : !! SOURCE
1789 :
1790 : subroutine gaussfit_param2_findsign()
1791 :
1792 : !Arguments -------------------------------
1793 : !Local variables-------------------------------
1794 : integer::ix,minx
1795 : real(dp)::dist,mindist,xx,yy
1796 :
1797 : ! *********************************************************************
1798 :
1799 : mindist=rpaw
1800 : do ix=1,nx
1801 : xx=x(ix)
1802 : dist=abs(raux-xx)
1803 : if(dist<mindist) then
1804 : mindist=dist
1805 : minx=ix
1806 : end if
1807 : end do
1808 : yy=y(minx)
1809 : sig=yy
1810 :
1811 : end subroutine gaussfit_param2_findsign
1812 :
1813 : end subroutine gaussfit_set_param2
1814 : !!***
1815 :
1816 : !----------------------------------------------------------------------
1817 :
1818 : !!****f* m_paw_gaussfit/gaussfit_set_param3
1819 : !! NAME
1820 : !! gaussfit_set_param3
1821 : !!
1822 : !! FUNCTION
1823 : !! Sets parameters for LSF
1824 : !!
1825 : !! INPUTS
1826 : !!
1827 : !! OUTPUT
1828 : !!
1829 : !! SOURCE
1830 :
1831 : subroutine gaussfit_set_param3(nterm,nparam,param,sep)
1832 :
1833 : !Arguments -------------------------------
1834 : integer,intent(in)::nterm,nparam
1835 : real(dp),intent(in)::sep
1836 : !real(dp),intent(in)::x(nx),y(nx)
1837 : real(dp),intent(out)::param(nparam)
1838 :
1839 : !Local variables-------------------------------
1840 : integer::ii,jj
1841 :
1842 : ! *********************************************************************
1843 :
1844 : param(:)=1.0d0
1845 : !
1846 : do jj=1,nterm
1847 : ii=jj+nterm
1848 : param(ii)=sep**(jj)
1849 : ! param(ii)=raux*real(i,dp)
1850 : end do
1851 : !
1852 : do jj=1,nterm
1853 : ii=jj+3*nterm
1854 : param(ii)=sep**(jj)
1855 : ! param(ii)=raux*real(i,dp)
1856 : end do
1857 :
1858 : end subroutine gaussfit_set_param3
1859 : !!***
1860 :
1861 : !----------------------------------------------------------------------
1862 :
1863 : !!****f* m_paw_gaussfit/gaussfit_set_param4
1864 : !! NAME
1865 : !! gaussfit_set_param4
1866 : !!
1867 : !! FUNCTION
1868 : !! Sets parameters for LSF
1869 : !!
1870 : !! INPUTS
1871 : !!
1872 : !! OUTPUT
1873 : !!
1874 : !! SOURCE
1875 :
1876 0 : subroutine gaussfit_set_param4(nparam,param)
1877 :
1878 : !Arguments -------------------------------
1879 : integer,intent(in)::nparam
1880 : real(dp),intent(out)::param(nparam)
1881 :
1882 : !Local variables-------------------------------
1883 :
1884 : ! *********************************************************************
1885 :
1886 0 : param(:)=1.0d0
1887 :
1888 : end subroutine gaussfit_set_param4
1889 : !!***
1890 :
1891 : !----------------------------------------------------------------------
1892 :
1893 : !!****f* m_paw_gaussfit/gaussfit_set_param5
1894 : !! NAME
1895 : !! gaussfit_set_param5
1896 : !!
1897 : !! FUNCTION
1898 : !! Sets parameters for LSF
1899 : !!
1900 : !! INPUTS
1901 : !!
1902 : !! OUTPUT
1903 : !!
1904 : !! SOURCE
1905 :
1906 0 : subroutine gaussfit_set_param5(nterm,nparam,nx,param,rpaw,y)
1907 :
1908 : !Arguments -------------------------------
1909 : integer,intent(in)::nterm,nparam,nx
1910 : real(dp),intent(in)::rpaw
1911 : real(dp),intent(in)::y(nx)
1912 : real(dp),intent(out)::param(nparam)
1913 :
1914 : !Local variables-------------------------------
1915 : integer::ix
1916 : real(dp)::raux,a1,r_c,m
1917 :
1918 : ! *********************************************************************
1919 :
1920 0 : param(:)=1.0d0
1921 : !
1922 : !alpha1
1923 : !a1=maxval( y(:),nx )
1924 : !maxval gives problems in some architectures:
1925 : a1=-9999999
1926 0 : do ix=1,nx
1927 0 : if(a1<y(ix)) a1=y(ix)
1928 : end do
1929 0 : param(1:nterm)=a1
1930 : !
1931 : !alpha2
1932 : !
1933 0 : r_c=rpaw+0.5d0 !paw sphere + a bit more.
1934 : !this is not arbitrary since it is an initial guess
1935 0 : m=0.01d0
1936 0 : raux=log(a1/m)/r_c**2
1937 0 : param(nterm+1:nterm*2)=raux
1938 :
1939 0 : end subroutine gaussfit_set_param5
1940 : !!***
1941 :
1942 : !----------------------------------------------------------------------
1943 :
1944 : !!****f* m_paw_gaussfit/gaussfit_constrains_init
1945 : !! NAME
1946 : !! gaussfit_constrains_init
1947 : !!
1948 : !! FUNCTION
1949 : !! Initialise constrains for LSF.
1950 : !! It will constrain the Gaussians width
1951 : !! It will also constraint the Delta use in the LSF algorithm, to jump slowly at each step.
1952 : !!
1953 : !! INPUTS
1954 : !!
1955 : !! OUTPUT
1956 : !!
1957 : !! SOURCE
1958 :
1959 0 : subroutine gaussfit_constrains_init(cons1,cons2,limit,nparam,nterm,nx,option,rpaw,y)
1960 :
1961 : !Arguments -------------------------------
1962 : integer,intent(in)::nterm,option,nparam,nx
1963 : integer,intent(out)::cons2(nparam)
1964 : real(dp),intent(in)::rpaw
1965 : real(dp),intent(in)::y(nx)
1966 : real(dp),intent(out)::cons1(nparam),limit(nparam)
1967 :
1968 : !Local variables-------------------------------
1969 : integer :: ix
1970 : real(dp)::rc,a1,mm,raux
1971 :
1972 : ! *********************************************************************
1973 :
1974 : !
1975 : !DEFAULT: no weight
1976 0 : cons1(:)=1.d0
1977 0 : limit(:)=0.d0
1978 0 : cons2=1
1979 : !
1980 0 : if(option==4) then
1981 0 : cons1(1:nterm)=0.2d0
1982 0 : cons1(nterm+1:nterm*2)=0.3d0
1983 : end if
1984 : !
1985 0 : if(option==4) then
1986 : ! parameters
1987 :
1988 : ! a1=maxval( y(:),nx)/real(nterm,dp)
1989 : ! maxval gives problems in some architectures:
1990 : a1=-9999999
1991 0 : do ix=1,nx
1992 0 : if(a1<y(ix)) a1=y(ix)
1993 : end do
1994 0 : a1=a1/real(nterm,dp)
1995 :
1996 0 : mm=0.01
1997 : !
1998 0 : rc=1.7d0*rpaw
1999 0 : raux=log(a1/mm)/rc**2
2000 : !
2001 : ! Constraint exponential of gaussians to be positive (multiplied by -1),
2002 : ! so that it decays to zero.
2003 : ! Constraint as well its value, so that it does not get too big
2004 : ! and it decays soon,
2005 : ! This prevents that a gaussian grows at a very large x value.
2006 0 : cons2(nterm+1:nterm*2)=restricted_and_positive
2007 0 : limit(nterm+1:nterm*2)=raux
2008 : end if
2009 :
2010 0 : end subroutine gaussfit_constrains_init
2011 : !!***
2012 :
2013 : !----------------------------------------------------------------------
2014 :
2015 : !!****f* m_paw_gaussfit/gaussfit_apply_constrains
2016 : !! NAME
2017 : !! gaussfit_apply_constrains
2018 : !!
2019 : !! FUNCTION
2020 : !! Apply constrains to get new set of parameters
2021 : !!
2022 : !! INPUTS
2023 : !!
2024 : !! OUTPUT
2025 : !!
2026 : !! SOURCE
2027 :
2028 0 : subroutine gaussfit_apply_constrains(const,limit,nparam,ioparams)
2029 :
2030 : !Arguments -------------------------------
2031 : integer,intent(in):: nparam
2032 : integer,intent(in):: const(nparam)
2033 : real(dp),intent(in):: limit(nparam)
2034 : real(dp),intent(inout):: ioparams(nparam)
2035 :
2036 : !Local variables-------------------------------
2037 : integer::ii
2038 :
2039 : ! *********************************************************************
2040 :
2041 0 : do ii=1,nparam
2042 0 : if(const(ii)==restricted .or. const(ii)==restricted_and_positive) then
2043 0 : if(ioparams(ii)<limit(ii)) ioparams(ii)=limit(ii)
2044 : end if
2045 0 : if(const(ii)==positive .or. const(ii)==restricted_and_positive ) ioparams(ii)=abs(ioparams(ii))
2046 :
2047 : end do
2048 :
2049 0 : end subroutine gaussfit_apply_constrains
2050 : !!***
2051 :
2052 : !----------------------------------------------------------------------
2053 :
2054 : !!****f* m_paw_gaussfit/gaussfit_projector
2055 : !! NAME
2056 : !! gaussfit_projector
2057 : !!
2058 : !! FUNCTION
2059 : !! Fit tproj to Gaussians
2060 : !!
2061 : !! INPUTS
2062 : !! basis_size= size of the PAW basis
2063 : !! orbitals= indicates the l quantum number for all orbitals.
2064 : !! rpaw= PAW radius
2065 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data
2066 : !! tproj= projectors
2067 : !! maxterm= maximum number of terms used to fit the projectors.
2068 : !! mparam= maximum number of parameters (Gaussian coefficients and factors) used.
2069 : !!
2070 : !! OUTPUT
2071 : !! nparam_array= number of parameters found.
2072 : !! param = parameters found (Gaussian coefficients and factors).
2073 : !!
2074 : !! NOTES
2075 : !! chisq=accuracy_p= sum_x abs(f(x)-y(x))/nx.
2076 : !! nx is the number of points, f(x) and y(x) are the fitted and original functions.
2077 : !!
2078 : !! SOURCE
2079 :
2080 0 : subroutine gaussfit_projector(basis_size,mparam,nparam_array,nterm_bounds,orbitals,param,pawrad,&
2081 0 : & rpaw,tproj,comm_mpi)
2082 :
2083 : !Arguments ------------------------------------
2084 : integer,intent(in) :: basis_size
2085 : integer,intent(in) :: orbitals(basis_size)
2086 : integer, optional,intent(in) :: comm_mpi
2087 : real(dp),intent(in) :: rpaw
2088 : type(pawrad_type),intent(in) :: pawrad
2089 : real(dp),intent(in) :: tproj(:,:)
2090 : integer,intent(in) :: mparam,nterm_bounds(2)
2091 : integer,intent(out) :: nparam_array(basis_size)
2092 : real(dp),intent(out) :: param(mparam,basis_size)
2093 0 : type(pawrad_type)::mesh_tmp
2094 :
2095 : !Local variables ------------------------------
2096 : integer :: ibasis,ierr,il,ir
2097 : integer :: msz1,msz2,option
2098 : real(dp) :: raux(1),rr(1)
2099 0 : real(dp),allocatable :: d2(:),tproj_tmp1(:),tproj_tmp2(:)
2100 : character(len=500) :: msg
2101 : character(80) :: outfile
2102 : !debug: uncomment
2103 : !integer::i,nterm ,unitp
2104 : !real(dp),allocatable::y(:)
2105 : !end debug
2106 :
2107 : !************************************************************************
2108 :
2109 0 : if(size(tproj,2)<basis_size) then
2110 0 : msg = 'wrong size for tproj in gaussfit_projector!'
2111 0 : LIBPAW_BUG(msg)
2112 : end if
2113 :
2114 0 : option=4 !see gaussfit_main
2115 0 : nparam_array(:)=0
2116 0 : msz1=min(pawrad_ifromr(pawrad,rpaw)+2,size(tproj,1))
2117 : !msz1=pawrad%mesh_size
2118 :
2119 : !Augment the mesh size
2120 : !this is to make the Gaussians go to zero after paw_radius
2121 : !This is done by creating a new pawrad objet: mesh_tmp
2122 : !Change to a linear grid:
2123 0 : mesh_tmp%mesh_type=1 !linear grid
2124 0 : mesh_tmp%rstep=0.0005 !very fine grid
2125 0 : msz2=ceiling(pawrad%rmax*two/mesh_tmp%rstep)
2126 : mesh_tmp%lstep=zero !only needed for log grids
2127 : call pawrad_init(mesh_tmp,mesh_size=msz2,mesh_type=mesh_tmp%mesh_type,&
2128 0 : & rstep=mesh_tmp%rstep,lstep=mesh_tmp%lstep)
2129 :
2130 0 : LIBPAW_ALLOCATE(tproj_tmp1,(msz1))
2131 0 : LIBPAW_ALLOCATE(d2,(msz1))
2132 0 : LIBPAW_ALLOCATE(tproj_tmp2,(msz2))
2133 :
2134 0 : do ibasis=1,basis_size
2135 :
2136 0 : write(msg,'(a," - Fitting wfn ",i4," to Gaussians")')ch10,ibasis
2137 0 : call wrtout(std_out, msg,'COLL')
2138 :
2139 0 : tproj_tmp1=zero; d2=zero; tproj_tmp2=zero
2140 :
2141 : ! take out r^il factor:
2142 : ! il=psps%indlmn(1,ilmn,itypat)
2143 0 : il=orbitals(ibasis)
2144 :
2145 0 : tproj_tmp1(2:msz1)=tproj(2:msz1,ibasis)/((pawrad%rad(2:msz1)+tol8)**(il))
2146 :
2147 : ! take out 1/r factor from eq.(3) of M. Torrent CMS 42, 337 (2008)
2148 : ! since: <phi|proj>=1 from atompaw, and phi=phi*r, alors proj=proj/r
2149 :
2150 0 : tproj_tmp1(2:msz1)=tproj_tmp1(2:msz1)/(pawrad%rad(2:msz1))
2151 0 : call pawrad_deducer0(tproj_tmp1(1:msz1),msz1,pawrad)
2152 :
2153 : ! splint to a different mesh:
2154 : ! get second derivative of tproj and store it
2155 : call paw_spline(pawrad%rad,tproj_tmp1(:),msz1,&
2156 0 : & zero,zero,d2)
2157 :
2158 0 : do ir=2,msz2
2159 0 : rr=mesh_tmp%rad(ir)
2160 0 : if( rr(1)-rpaw > tol8 ) then
2161 : !after rpaw projectors are zero
2162 0 : raux=zero
2163 : else
2164 : call paw_splint(msz1,pawrad%rad,&
2165 : & tproj_tmp1(:),d2(:),&
2166 0 : & 1,rr,raux,ierr=ierr)
2167 : end if
2168 0 : tproj_tmp2(ir)=raux(1)
2169 : end do
2170 :
2171 : ! Obtain the name for the output file
2172 0 : if(ibasis<10) then
2173 0 : write(outfile,'("wfn",i1,".fit")')ibasis
2174 0 : elseif(ibasis<100) then
2175 0 : write(outfile,'("wfn",i2,".fit")')ibasis
2176 0 : write(msg,'(a,a,a,a)')ch10,&
2177 0 : & "ib (basis index) is too big!",ch10,&
2178 0 : & "Action: check your pseudopotentials"
2179 0 : LIBPAW_BUG(msg)
2180 : end if
2181 :
2182 0 : if(present(comm_mpi)) then
2183 : call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,&
2184 0 : & param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2,comm_mpi)
2185 : else
2186 : call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,&
2187 0 : & param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2)
2188 : end if
2189 :
2190 : ! check
2191 : ! LIBPAW_ALLOCATE(y,(mesh_tmp%mesh_size))
2192 : ! nterm=nparam_array(ibasis)/4
2193 : ! call calcgaussc4(nparam_array(ibasis),nterm,mesh_tmp%mesh_size,1,param(:,ibasis),&
2194 : ! & mesh_tmp%rad,y)
2195 : ! ! unitp=600+ibasis
2196 : ! ! do ir=1,mesh_tmp%mesh_size
2197 : ! ! write(unitp,'(2(f16.7,x,f16.7))')mesh_tmp%rad(ir),y(ir)
2198 : ! ! end do
2199 : ! LIBPAW_DEALLOCATE(y)
2200 :
2201 : end do
2202 :
2203 : !Deallocate
2204 0 : call pawrad_free(mesh_tmp)
2205 0 : LIBPAW_DEALLOCATE(tproj_tmp1)
2206 0 : LIBPAW_DEALLOCATE(tproj_tmp2)
2207 0 : LIBPAW_DEALLOCATE(d2)
2208 :
2209 0 : end subroutine gaussfit_projector
2210 : !!***
2211 :
2212 : !----------------------------------------------------------------------
2213 :
2214 : end module m_paw_gaussfit
2215 : !!***
|
