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