Branch data Line data Source code
1 : : !!****m* ABINIT/m_rec
2 : : !! NAME
3 : : !! m_rec
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides some functions applied to the
7 : : !! recursion structured datatype recursion_type.
8 : : !! It includes also some function used to change some variables
9 : : !! of recursion_type
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2002-2026 ABINIT group (MMancini)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
17 : : !!
18 : : !! NOTES
19 : : !!
20 : : !! * Routines tagged with "@type_name" are strongly connected to the definition of the data type.
21 : : !! Strongly connected means that the proper functioning of the implementation relies on the
22 : : !! assumption that the tagged procedure is consistent with the type declaration.
23 : : !! Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
24 : : !! that all the strongly connected routines are changed accordingly to accomodate the modification of the data type.
25 : : !! Typical examples of strongly connected routines are creation, destruction or reset methods.
26 : : !!
27 : : !!
28 : : !! SOURCE
29 : :
30 : : #if defined HAVE_CONFIG_H
31 : : #include "config.h"
32 : : #endif
33 : :
34 : : #include "abi_common.h"
35 : :
36 : : module m_rec
37 : :
38 : : use defs_basis
39 : : use defs_rectypes
40 : : use m_abicore
41 : : use m_errors
42 : : use m_xmpi
43 : : use m_sort
44 : : use m_dtset
45 : :
46 : : use defs_datatypes, only : pseudopotential_type
47 : : use defs_abitypes, only : mpi_type
48 : : use m_exp_mat, only : exp_mat
49 : : use m_numeric_tools, only : set2unit
50 : : use m_special_funcs, only : gamma_function
51 : : use m_pawfgr, only : pawfgr_nullify, indgrid, pawfgr_destroy
52 : : use m_paw_sphharm, only : initylmr
53 : : use m_time, only : timab
54 : : use m_rec_tools, only : get_pt0_pt1
55 : : use m_per_cond, only : per_cond
56 : : #ifdef HAVE_GPU_CUDA
57 : : use m_hidecudarec, only : InitRecGPU, CleanRecGPU
58 : : #endif
59 : :
60 : : implicit none
61 : :
62 : : private :: &
63 : : find_maxmin_proc, & !--To calculate max and min pt for any cpu
64 : : H_D_distrib
65 : :
66 : : public :: &
67 : : InitRec, & !--Main creation method.
68 : : Init_MetricRec, & !--To Initalize the inf. metric in recursion
69 : : Init_nlpspRec, & !--Main creation method for non-local part.
70 : : CleanRec, & !--deallocate all pointers.
71 : : Calcnrec, & !--calculates the new min_nrec
72 : : cpu_distribution !--Regulates the work load on cpu-gpu
73 : : CONTAINS !===========================================================
74 : : !!***
75 : :
76 : : !!****f* m_rec/H_D_distrib
77 : : !! NAME
78 : : !! H_D_distrib
79 : : !!
80 : : !! FUNCTION
81 : : !! Calculate the number of point,GPU,for any proc
82 : : !!
83 : : !! INPUTS
84 : : !! rset<recursion_type>= recursion variables
85 : : !! cpu (-1 if there are not gpu)
86 : : !! nfft=nuber of point of the fine grid
87 : : !! ngfftrec=nuber of point of one edge of the coarse grid
88 : : !! gratio=recgratio ratio between the fine and coarse grid
89 : : !! beta_coeff=estimated time ratio between CPU_time and GPU_time
90 : : !!
91 : : !! OUTPUT
92 : : !! proc_pt_dev(2,0:nproc-1) which device and how many points
93 : : !! that proc has to compute: proc_pt_dev(1,iproc) which device
94 : : !! associated to proc i (-1 if none), proc_pt_dev(2,iproc) how
95 : : !! many points
96 : : !!
97 : : !! SOURCE
98 : :
99 : 6 : subroutine H_D_distrib(rset,nfft,gratio,proc_pt_dev,beta_coeff)
100 : :
101 : : !Arguments ------------------------------------
102 : : integer, intent(in) :: nfft,gratio
103 : : real(dp),intent(in) :: beta_coeff
104 : : integer,pointer :: proc_pt_dev(:,:)
105 : : type(recursion_type),intent(inout) :: rset
106 : : !Local ---------------------------
107 : : integer :: me,icpu,resto,ntot,ngpu
108 : : integer :: n_per_cpu,n_per_gpu
109 : : character(500) :: msg
110 : : #ifdef HAVE_GPU_CUDA
111 : : integer,pointer :: ndev(:)
112 : : #else
113 : 12 : integer :: ndev(0:rset%mpi%nproc-1)
114 : : #endif
115 : : ! *********************************************************************
116 : :
117 : :
118 : : #ifdef HAVE_GPU_CUDA
119 : : ndev => rset%GPU%map
120 : : #else
121 [ + + ]: 24 : ndev = -1
122 : : #endif
123 : :
124 : 6 : me = rset%mpi%me
125 : 6 : ntot = nfft/(gratio*gratio*gratio)
126 : 6 : ngpu = rset%ngpu
127 : :
128 : : !--If sequential code all points are computed by the proc 0
129 [ + + ]: 6 : if(rset%mpi%nproc ==1) then
130 : 2 : proc_pt_dev(1,0) = ndev(0)
131 : 2 : proc_pt_dev(2,0) = ntot
132 : 2 : return
133 : : end if
134 : :
135 : : !--Number of points for any cpu
136 : 4 : n_per_cpu = int(int(ntot/(rset%mpi%nproc+ngpu*(beta_coeff-1.d0))))
137 : 4 : n_per_gpu = int(n_per_cpu*beta_coeff)
138 : : !write(std_out,*)'n_per_cpu',n_per_cpu
139 : : !write(std_out,*)'rset%GPU%map',rset%GPU%map
140 [ + + ]: 20 : do icpu=0,rset%mpi%nproc-1
141 : 16 : proc_pt_dev(1,icpu) = ndev(icpu)
142 : 16 : proc_pt_dev(2,icpu) = n_per_cpu
143 [ - + ]: 20 : if(ndev(icpu)>-1) proc_pt_dev(2,icpu) = n_per_gpu
144 : : end do
145 : :
146 : : !--Distribute the rest
147 [ + + ]: 20 : resto = ntot-sum(proc_pt_dev(2,:))
148 : 4 : icpu = 0
149 : : !write(std_out,*)'rest',resto,ngpu
150 [ - + ]: 4 : if(resto>0) then
151 [ # # ]: 0 : if(ngpu/=0) then
152 : : !--distribute rest only on GPU
153 [ # # ]: 0 : do while(resto/=0)
154 [ # # ]: 0 : if(proc_pt_dev(1,icpu)>-1) then
155 : 0 : proc_pt_dev(2,icpu) = proc_pt_dev(2,icpu)+1
156 : 0 : resto = resto-1
157 : : endif
158 : 0 : icpu = mod(icpu+1,rset%mpi%nproc)
159 : : enddo
160 : : else
161 : : !--distribute rest on all CPU
162 [ # # ]: 0 : do while(resto/=0)
163 : 0 : proc_pt_dev(2,icpu) = proc_pt_dev(2,icpu)+1
164 : 0 : resto = resto-1
165 : 0 : icpu = mod(icpu+1,rset%mpi%nproc)
166 : : enddo
167 : : return
168 : : endif
169 : : endif
170 : :
171 : : !--Printing GPU and load distribution on procs
172 : : write(msg,'(3a)')&
173 : 4 : & ' -Load on procs------------',ch10,&
174 : 8 : & ' me device points'
175 : 4 : call wrtout(std_out,msg,'COLL')
176 [ + + ]: 20 : do icpu=0,rset%mpi%nproc-1
177 : 16 : write(msg,'(i5,i8,i14)') icpu,proc_pt_dev(:,icpu);
178 : 20 : call wrtout(std_out,msg,'COLL')
179 : : end do
180 : :
181 : : end subroutine H_D_distrib
182 : : !!***
183 : :
184 : :
185 : :
186 : : !!****f* m_rec/find_maxmin_proc
187 : : !! NAME
188 : : !! find_maxmin_proc
189 : : !!
190 : : !! FUNCTION
191 : : !! To calculate max and min pt for any cpu, it is useful for
192 : : !! recgratio!=1
193 : : !!
194 : : !! INPUTS
195 : : !! nproc = number of procs
196 : : !! me = identity of the proc
197 : : !! ngfft(3) = fine grid (corresponds to dtset%ngfft(1:3))
198 : : !! proc_pt_dev(2,0:nproc-1) which device and how many points
199 : : !! recpar%npt = number of points computed by the proc me (see side effects)
200 : : !!
201 : : !! OUTPUT
202 : : !! recpar%pt0<type(vec_int)>=Intial point for this proc in x,y,z
203 : : !! recpar%pt1<type(vec_int)>=Final point for this proc in x,y,z
204 : : !! recpar%min_pt=Intial point for this proc
205 : : !! recpar%max_pt=Final point for this proc
206 : : !!
207 : : !! SIDE EFFECTS
208 : : !! recpar%ntranche=number of pts computed by the proc me on the fine grid.
209 : : !!
210 : : !!
211 : : !! So when recgratio!=1, ntranche will not correspond to the npt!
212 : : !!
213 : : !! SOURCE
214 : :
215 : 6 : subroutine find_maxmin_proc(recpar,nproc,me,gratio,ngfft,proc_pt_dev)
216 : :
217 : : !Arguments ------------------------------------
218 : : integer,intent(in) :: nproc,me,gratio
219 : : integer,intent(in) :: ngfft(3)
220 : : type(recparall_type),intent(inout) :: recpar
221 : : integer,pointer :: proc_pt_dev(:,:)
222 : : !Local ---------------------------
223 : : integer :: pointoncpu
224 : : integer :: nfft,ntot,ii
225 : : integer :: inf,sup
226 : 12 : integer :: proc_limit(0:nproc-1)
227 : : ! *********************************************************************
228 : : ! write(std_out,*)'start find_maxmin_proc'
229 : 6 : recpar%npt = 0
230 : : nfft = product(ngfft)
231 : 6 : ntot = nfft/(gratio*gratio*gratio)
232 : 6 : pointoncpu = ntot/nproc
233 : :
234 [ + - + + : 90 : proc_limit = (/(sum(proc_pt_dev(2,:ii)),ii=0,nproc-1)/)
+ + + + ]
235 : :
236 [ + + ]: 6 : if(gratio==1)then
237 : 5 : recpar%ntranche = proc_limit(me)
238 [ + + ]: 5 : if(me/=0) recpar%ntranche = recpar%ntranche-proc_limit(me-1)
239 : : endif
240 : :
241 : 6 : inf=0
242 [ + + ]: 6 : if(me/=0) inf = proc_limit(me-1)
243 : 6 : sup = proc_limit(me)
244 : :
245 : :
246 : 6 : call get_pt0_pt1(ngfft,gratio,inf,sup,recpar)
247 : :
248 : 6 : recpar%npt = sup-inf
249 : :
250 : : !write(std_out,*)'exit find_maxmin_proc'
251 : 6 : end subroutine find_maxmin_proc
252 : : !!***
253 : :
254 : : !!****f* m_rec/cpu_distribution
255 : : !! NAME
256 : : !! cpu_distribution
257 : : !!
258 : : !! FUNCTION
259 : : !! Calculate the number of point,GPU,for any proc
260 : : !!
261 : : !! INPUTS
262 : : !! ngfft(3)=nuber of point of the grid
263 : : !! gratio=recgratio ratio between the fine and coarse grid
264 : : !! beta_coeff=estimated time ratio between CPU_time and GPU_time
265 : : !! calc_type=if 0 takes the possible max for nptrec (to test the
266 : : !! completly full graphic card). 1 after test to calculate the min
267 : : !! possible value for nptrec
268 : : !!
269 : : !! OUTPUT
270 : : !!
271 : : !! SOURCE
272 : :
273 : 6 : subroutine cpu_distribution(gratio,rset,ngfft,beta_coeff,calc_type)
274 : :
275 : : !Arguments ------------------------------------
276 : : integer,intent(in) :: gratio,calc_type
277 : : real(dp),intent(in) :: beta_coeff
278 : : integer,intent(in) :: ngfft(3)
279 : : type(recursion_type),intent(inout),target :: rset
280 : : !Local ---------------------------
281 : : integer :: ii,nfft,ierr
282 : 6 : integer,pointer :: proc_pt_dev(:,:)
283 : : type(recparall_type),pointer :: recpar
284 : : character(500) :: msg
285 : : ! *********************************************************************
286 : :
287 : : ! write(std_out,*)'start cpu_distribution'
288 : :
289 : : nullify(proc_pt_dev)
290 [ - + + - : 18 : ABI_MALLOC(proc_pt_dev,(2,0:rset%mpi%nproc-1))
+ - + - -
+ - + ]
291 : :
292 [ + + ]: 24 : nfft = product(ngfft)
293 : 6 : call H_D_distrib(rset,nfft,gratio,proc_pt_dev,beta_coeff)
294 : :
295 : 6 : nullify(recpar)
296 [ + - ]: 6 : if(rset%load == 0)then
297 [ - + + - : 18 : ABI_MALLOC(rset%par%displs,(0:rset%mpi%nproc-1))
+ - + - -
+ - + -
+ ]
298 [ - + + - : 18 : ABI_MALLOC(rset%par%vcount,(0:rset%mpi%nproc-1))
+ - + - -
+ - + -
+ ]
299 : 6 : recpar => rset%par
300 : : #if defined HAVE_GPU_CUDA
301 : : else
302 : : if(rset%tp==4)then
303 : : if(.not. allocated(rset%GPU%par%displs)) then
304 : : ABI_MALLOC(rset%GPU%par%displs,(0:rset%mpi%nproc-1))
305 : : ABI_MALLOC(rset%GPU%par%vcount,(0:rset%mpi%nproc-1))
306 : : end if
307 : : endif
308 : : recpar => rset%GPU%par
309 : : #endif
310 : : endif
311 : :
312 : 6 : recpar%ntranche = nfft/(rset%mpi%nproc)!equipartitioned point
313 : :
314 : : call find_maxmin_proc(recpar,rset%mpi%nproc,&
315 : 6 : & rset%mpi%me,gratio,ngfft,proc_pt_dev)
316 : :
317 [ + + ]: 24 : recpar%vcount = 0
318 [ + - ]: 6 : if(rset%load==0)then
319 : 6 : recpar%vcount(rset%mpi%me) = recpar%ntranche
320 : : else
321 : 0 : recpar%vcount(rset%mpi%me) = recpar%npt
322 : : endif
323 : :
324 : 6 : call xmpi_sum(recpar%vcount,rset%mpi%comm_bandfft,ierr)
325 : :
326 [ + + ]: 24 : recpar%displs = 0
327 [ + + + - : 58 : if(rset%mpi%nproc>1) recpar%displs(1:) = (/(sum(recpar%vcount(:ii)),ii=0,rset%mpi%nproc-2)/)
+ + + + +
+ ]
328 : :
329 : : !--INITALIZATION OF CUDA FOR RECURSION
330 : : #if defined HAVE_GPU_CUDA
331 : : if(rset%load == 0) rset%GPU%par = rset%par
332 : : call InitRecGPU(rset,nfft,gratio,rset%GPU%map(rset%mpi%me),calc_type)
333 : : #else
334 : : ierr = calc_type !only of abirule when there is not HAVE_GPU_CUDA
335 : : #endif
336 : :
337 : :
338 : : ! if(rset%debug ) then
339 : : write(msg,'(a,i7,2(2a,3i7),8(2a,i7),2(2a,3i7),(2a,e14.6))')&
340 : 6 : & ' me ', rset%mpi%me,ch10,&
341 : 6 : & ' ngfft ', ngfft(1:3),ch10,&
342 : 6 : & ' ngfftrec ', rset%ngfftrec(1:3),ch10,&
343 : 6 : & ' load ', rset%load,ch10,&
344 : 6 : & ' ntranche ', recpar%ntranche,ch10,&
345 : 6 : & ' min_pt ', recpar%min_pt,ch10,&
346 : 6 : & ' max_pt ', recpar%max_pt,ch10,&
347 : 6 : & ' rset%mpi%nproc ', rset%mpi%nproc,ch10,&
348 : 6 : & ' rset%mpi%nproc_fft ', rset%mpi%nproc_fft,ch10,&
349 : 6 : & ' dtset%ngfft(10) ', rset%ngfftrec(10),ch10,&
350 : 6 : & ' recpar%npt ', recpar%npt,ch10,&
351 : 6 : & ' recpar%pt0 ', recpar%pt0%x,recpar%pt0%y,recpar%pt0%z,ch10,&
352 : 6 : & ' recpar%pt1 ', recpar%pt1%x,recpar%pt1%y,recpar%pt1%z,ch10,&
353 : 12 : & ' grid step ', rset%inf%tr(1)
354 : 6 : call wrtout(std_out,msg,'PERS')
355 : : #if defined HAVE_GPU_CUDA
356 : : write(msg,'(a,i7,2(2a,i7),a)')&
357 : : & ' rset%ngp ', rset%ngpu,ch10,&
358 : : & ' gpudevice ', rset%gpudevice,ch10,&
359 : : & ' nptrec ', rset%GPU%nptrec,ch10
360 : : call wrtout(std_out,msg,'PERS')
361 : : #endif
362 : : ! write(std_out,*)'display',recpar%displs
363 : : ! write(std_out,*)'vcount',recpar%vcount
364 : : ! end if
365 : :
366 : :
367 : 6 : nullify(recpar)
368 [ + - ]: 6 : if(associated(proc_pt_dev)) then
369 : 6 : ABI_FREE(proc_pt_dev)
370 : : end if
371 : :
372 : : ! write(std_out,*)'exit from cpu_distribution'
373 : 6 : end subroutine cpu_distribution
374 : : !!***
375 : :
376 : :
377 : : !!****f* m_rec/InitRec
378 : : !! NAME
379 : : !! InitRec
380 : : !!
381 : : !! FUNCTION
382 : : !! Initialise the rset<recursion_type>=Data type concerning recursion.
383 : : !!
384 : : !! INPUTS
385 : : !! dtset <type(dataset_type)>=all input variables in this dataset
386 : : !! mpi_ab <type(mpi_type)=MPI-parallelisation information
387 : : !! mproj=0 if psp is only local
388 : : !!
389 : : !! SIDE EFFECTS
390 : : !! All pointers set to null().
391 : : !!
392 : : !!
393 : : !! SOURCE
394 : :
395 : 6 : subroutine InitRec(dtset,mpi_ab,rset,rmet,mproj)
396 : :
397 : : #ifdef HAVE_GPU_CUDA
398 : : use m_gpu_detect,only :get_topo,find_set_gpu
399 : : use m_hidecudarec,only :InitRecGPU_0
400 : : #include "cuda_common.h"
401 : : #endif
402 : :
403 : : !Arguments ------------------------------------
404 : : ! scalars
405 : : integer,intent(in) :: mproj
406 : : type(dataset_type),intent(in) :: dtset
407 : : type(MPI_type),intent(in),target :: mpi_ab
408 : : type(recursion_type),intent(inout) :: rset
409 : : real(dp),intent(in) :: rmet(3,3)
410 : : ! arrays
411 : : !Local ---------------------------
412 : : integer :: ii
413 : : real(dp) :: beta,rtrotter
414 : : #if defined HAVE_GPU_CUDA
415 : : character(500) :: msg
416 : : #endif
417 : : ! *********************************************************************
418 : : ! @recursion_type
419 : : ! @pawfgr_type
420 : :
421 : : !--Initialisation
422 : 6 : beta = one/dtset%tsmear !--Inverse of temperature
423 : : !--Rewriting the trotter parameter
424 : 6 : rtrotter = max(half,real(dtset%recptrott,dp))
425 : :
426 : 6 : rset%debug= (dtset%prtvol==-7)
427 : 6 : rset%quitrec = 0
428 : 6 : rset%min_nrec = dtset%recnrec
429 : 6 : rset%efermi = dtset%recefermi !initial guess for fermie
430 : :
431 : 6 : rset%nfftrec = 0
432 [ + + ]: 114 : rset%ngfftrec = 0
433 : :
434 : 6 : rset%tronc = .False.
435 : :
436 : 6 : rset%mpi => mpi_ab
437 : :
438 : : !--Are all pseudo-potentials local?
439 : 6 : rset%nl%nlpsp = (mproj/=0)
440 : : !--Some initialisation concerning the metrics
441 : : ! If non-local psps then it allocates the atoms positions
442 : : ! on the grid
443 [ - + ]: 6 : if(rset%nl%nlpsp) then
444 [ # # # # : 0 : ABI_MALLOC(rset%inf%gcart,(3,dtset%natom))
# # # # #
# # # #
# ]
445 : : else
446 [ - + - + ]: 6 : ABI_MALLOC(rset%inf%gcart,(0,0))
447 : : end if
448 [ - + - - ]: 6 : rset%inf%gcart = 0
449 : :
450 : : !----------------------------------------------------------
451 : : !--TRONCATION OF THE BOX
452 : : !! determines new dimensions the method is similar to the one used
453 : : !! in getng (except that ecut and xboxcutmin give no constraint,
454 : : !! and symmetries are not handled)
455 : :
456 : 6 : call getngrec(dtset%ngfft,rmet,rset%ngfftrec,rset%nfftrec,dtset%recrcut,0.25d0*sqrt(beta/rtrotter),rset%tronc)
457 : : ! 1/4*sqrt(beta/trotter) for guess - should be modified
458 : :
459 : : !------------------------------------------------------------
460 : : !--DETERMINING WHICH POINT WILL COMPUTE THAT PROC
461 : : !----------------------------------------------------------
462 : : !--Paralelism using the band communicator (not used in the recursion)
463 : : !--Distribution on procs with cuda
464 : :
465 : :
466 : 6 : rset%ngpu = 0 !--Initial guess no GPU at all
467 : 6 : rset%gpudevice = -1 !--Inital guess no GPU associated
468 : 6 : rset%load = 0 !--Inital homogeneous work load
469 : 6 : rset%tp = 0 !--Initial guess 1 cpu, 0 gpu
470 : :
471 : :
472 : : #ifdef HAVE_GPU_CUDA
473 : : !--Initialise GPU variables for recursion
474 : : call InitRecGPU_0(rset%GPU,mpi_ab)
475 : :
476 : : !--Get the distribution of GPUs on CPUs
477 : : call find_set_gpu(mpi_ab%nproc,mpi_ab%comm_bandfft,rset%GPU%map,rset%ngpu)
478 : :
479 : : !--Get the topology of the machine
480 : : call get_topo(rset%mpi%nproc,rset%ngpu,rset%tp)
481 : : if(rset%tp>4)then
482 : : msg = 'm_rec: number of gpu>number of cpu is not implemented'
483 : : ABI_ERROR(msg)
484 : : endif
485 : : ! rset%tp = 0;if(rset%mpi%nproc>1)rset%tp = 1
486 : : ! rset%ngpu = 0; rset%GPU%map=-1
487 : : !--For the moment cuda doesnt take into account non-local psp
488 : : if(rset%nl%nlpsp) then
489 : : rset%tp = 0;if(rset%mpi%nproc>1)rset%tp = 1
490 : : rset%GPU%map = -1
491 : : endif
492 : : #else
493 [ + + ]: 6 : if(rset%mpi%nproc>1)rset%tp = 1
494 : : #endif
495 : :
496 : : !--Basic initialization for recursion metric (only needed for printing)
497 [ + + ]: 24 : do ii=1,3
498 [ + + ]: 78 : rset%inf%rmet(ii,:) = rmet(ii,:)/(real(dtset%ngfft(1:3)*dtset%ngfft(ii),dp))
499 : : end do
500 [ + + ]: 24 : rset%inf%tr(:) = sqrt((/(rset%inf%rmet(ii,ii),ii=1,3)/)) !grid step
501 : :
502 : : !--Compute the work loqd distribution on devices (gpu,cpu)
503 : 6 : call cpu_distribution(dtset%recgratio,rset,dtset%ngfft(:3),1.d0,0)
504 : :
505 : : !------------------------------------------------------------
506 : : !--DEFINITION VARIABLE COARSE-FINE GRID TO USE TRANSGRID--INGRID FUNCTIONS
507 : 6 : call pawfgr_nullify(rset%pawfgr)
508 : : !if coarse grid is used
509 [ + + ]: 6 : if (dtset%recgratio>1) then
510 : : !fine grid--
511 : 1 : rset%pawfgr%mgfft = 0
512 [ + + ]: 4 : rset%pawfgr%nfft = product(dtset%ngfft(1:3))
513 [ + + ]: 19 : rset%pawfgr%ngfft(:) = dtset%ngfft(:)
514 [ + + ]: 4 : rset%pawfgr%ngfft(9:11)=(/0,1,0/)
515 [ + + ]: 3 : rset%pawfgr%ngfft(12:13)= dtset%ngfft(2:3)
516 : : !coarse grid--
517 : 1 : rset%pawfgr%mgfftc = 0
518 [ + + ]: 19 : rset%pawfgr%ngfftc(:) = rset%pawfgr%ngfft(:)
519 [ + + - + ]: 4 : rset%pawfgr%ngfftc(:3) = floor(real(dtset%ngfft(:3)+1,dp)/real(dtset%recgratio,dp))
520 [ + + ]: 4 : rset%pawfgr%nfftc = product(rset%pawfgr%ngfftc(1:3))
521 : :
522 : 1 : rset%pawfgr%usefinegrid = 1
523 [ - + + - : 3 : ABI_MALLOC(rset%pawfgr%fintocoa,(rset%pawfgr%nfft))
+ - + - -
+ - + ]
524 [ - + + - : 3 : ABI_MALLOC(rset%pawfgr%coatofin,(rset%pawfgr%nfftc))
+ - + - -
+ - + ]
525 : : call indgrid(rset%pawfgr%coatofin,rset%pawfgr%fintocoa,&
526 : : rset%pawfgr%nfftc,rset%pawfgr%nfft,&
527 : 1 : rset%pawfgr%ngfftc,rset%pawfgr%ngfft)
528 : :
529 : : else
530 : 5 : rset%pawfgr%mgfft = 0
531 [ + + ]: 95 : rset%pawfgr%ngfft = 0
532 : 5 : rset%pawfgr%mgfftc = 0
533 : :
534 : 5 : rset%pawfgr%usefinegrid = 0
535 : : end if
536 : :
537 : :
538 : 6 : end subroutine InitRec
539 : : !!***
540 : :
541 : : !!****f* m_rec/Init_MetricRec
542 : : !! NAME
543 : : !! Init_MetricRec
544 : : !!
545 : : !! FUNCTION
546 : : !! Initialise the rset<recursion_type>=Data type concerning recursion.
547 : : !! In particular, the information on the infinitesimal metric.
548 : : !! Also other variable are initialized
549 : : !!
550 : : !! INPUTS
551 : : !! rmet: metrics
552 : : !! ucvol=unit cell volume in bohr**3.
553 : : !! ngfft(1:3)=fine grid used in recursion
554 : : !! rprimd=Real space PRIMitive translations, Dimensional
555 : : !! xred=vectors (X) of atom positions in reduced coordinates
556 : : !! natom=number of atoms
557 : : !! debug=debug variable
558 : : !!
559 : : !! OUTPUT
560 : : !! metrec <type(metricrec_type)>= infinitesimal metrics used in recursion
561 : : !!
562 : : !! SOURCE
563 : :
564 : 11 : subroutine Init_MetricRec(metrec,nlpsp,rmet,ucvol,rprimd,xred,ngfft,natom,debug)
565 : :
566 : : !Arguments ------------------------------------
567 : : !scalars
568 : : integer, intent(in) ::natom
569 : : real(dp), intent(in) :: ucvol
570 : : logical,intent(in) ::nlpsp,debug
571 : : type(metricrec_type),intent(inout) :: metrec
572 : : !arrays
573 : : integer,intent(in) :: ngfft(3)
574 : : real(dp),intent(in) :: rmet(3,3),rprimd(3,3),xred(3,natom)
575 : :
576 : : !Local ---------------------------
577 : : integer :: ii
578 : 22 : real(dp) :: xcart(3,natom)
579 : : character(500) :: msg
580 : : ! *********************************************************************
581 : :
582 : : !--Intialisation of variables concerning the infinitesimal metric
583 [ + + ]: 44 : do ii=1,3
584 [ + + ]: 143 : metrec%rmet(ii,:) = rmet(ii,:)/(real(ngfft(1:3)*ngfft(ii),dp))
585 : : end do
586 [ + + ]: 44 : metrec%ucvol = ucvol/real(product(ngfft(1:3)),dp)
587 [ + + ]: 44 : metrec%tr(:) = sqrt((/(metrec%rmet(ii,ii),ii=1,3)/)) !grid step
588 : :
589 : : !--Initialisation of others variables
590 : : !--In non-loc-psp case: calculate the position of ions and conversion factor
591 [ - + ]: 11 : if(nlpsp) then
592 [ # # ]: 0 : do ii = 1,natom
593 [ # # # # : 0 : xcart(:,ii) = matmul(rprimd(:,:),xred(:,ii))
# # ]
594 : : end do
595 [ # # # # : 0 : metrec%gcart(:,:) = per_cond(natom,xcart,ngfft(1:3),metrec%tr(:))
# # ]
596 [ # # ]: 0 : if(debug) then
597 [ # # ]: 0 : do ii=1,natom
598 : 0 : write (msg,'(a,3f8.2)')'xcart=',xcart(:,ii)
599 : 0 : call wrtout(std_out,msg,'COLL')
600 : 0 : write (msg,'(a,3i4)')'gcart=',metrec%gcart(:,ii)
601 : 0 : call wrtout(std_out,msg,'COLL')
602 : : end do
603 : : end if
604 : : end if
605 : :
606 : 11 : end subroutine Init_MetricRec
607 : : !!***
608 : :
609 : : !!****f* m_rec/Init_nlpspRec
610 : : !! NAME
611 : : !! Init_nlpspRec
612 : : !!
613 : : !! FUNCTION
614 : : !! Initialise the rset<recursion_type>=Data type concerning recursion.
615 : : !! In particular, the non-local part of pseudo-potential.
616 : : !!
617 : : !! INPUTS
618 : : !! tempe=temperature
619 : : !! psps <type(pseudopotential_type)>=variables related to pseudo-potentials
620 : : !! metrec <type(metricrec_type)>=infinitesimal metrics used in recursion
621 : : !! ngfftrec(18)=Number of Grid points for Fast Fourier Transform for
622 : : !! Recursion (truncated box, if different from ngfft)
623 : : !! debug=debug variable
624 : : !!
625 : : !! SIDE EFFECTS
626 : : !! nlrec <type(nlpsprec_type)>=pseudo-potentials information for recursion
627 : : !!
628 : : !! SOURCE
629 : :
630 : 6 : subroutine Init_nlpspRec(tempe,psps,nlrec,metrec,ngfftrec,debug)
631 : :
632 : : !Arguments ------------------------------------
633 : : ! scalars
634 : : logical,intent(in) :: debug
635 : : real(dp), intent(in) :: tempe
636 : : type(pseudopotential_type),intent(in) ::psps
637 : : type(metricrec_type),intent(inout) :: metrec
638 : : type(nlpsprec_type),intent(inout) :: nlrec
639 : : ! arrays
640 : : integer,intent(in) :: ngfftrec(18)
641 : : !Local ---------------------------
642 : : integer :: ii,jj
643 : : character(500) :: msg
644 : : ! *********************************************************************
645 : : !!--Routine for the calcul of the non-local pseudo
646 [ + + + - : 12 : if(any(psps%pspcod/=3) .and. nlrec%nlpsp ) then
- + - - ]
647 : 0 : msg = "The non-local part of psp is used in Recursion only for hgh-psp"
648 : 0 : ABI_WARNING(msg)
649 : 0 : nlrec%nlpsp = .False.
650 [ # # ]: 0 : if (allocated(metrec%gcart)) then
651 : 0 : ABI_FREE(metrec%gcart)
652 : : end if
653 : : end if
654 : :
655 [ + - - + : 6 : if(any(psps%pspcod==3) .and. nlrec%nlpsp) then
+ - - + ]
656 : :
657 : 0 : nlrec%nlpsp = .True.
658 : 0 : nlrec%npsp = psps%npsp
659 [ # # # # ]: 0 : nlrec%lmnmax = count(psps%indlmn(3,:,psps%npsp)/=0)
660 [ # # # # : 0 : ABI_MALLOC(nlrec%mat_exp_psp_nl,(3,3,psps%mpsang,psps%npsp))
# # # # #
# # # # #
# # # # ]
661 [ # # # # : 0 : ABI_MALLOC(nlrec%eival,(3,psps%mpsang,psps%npsp))
# # # # #
# # # # #
# # # # ]
662 [ # # # # : 0 : ABI_MALLOC(nlrec%eivec,(3,3,psps%mpsang,psps%npsp))
# # # # #
# # # #
# ]
663 [ # # # # : 0 : ABI_MALLOC(nlrec%pspinfo,(psps%mpsang,psps%npsp))
# # # # #
# # # # #
# # # # ]
664 [ # # # # : 0 : ABI_MALLOC(nlrec%radii,(psps%mpsang,psps%npsp))
# # # # #
# # # # #
# # # # ]
665 [ # # # # : 0 : ABI_MALLOC(nlrec%indlmn,(6,nlrec%lmnmax,psps%npsp))
# # # # #
# # # # #
# # # # ]
666 [ # # # # : 0 : nlrec%indlmn(:,:,:) = psps%indlmn(:,:nlrec%lmnmax,:)
# # ]
667 [ # # # # : 0 : nlrec%mat_exp_psp_nl(:,:,:,:) = zero
# # # # ]
668 [ # # # # : 0 : nlrec%eivec(:,:,:,:) = zero
# # # # ]
669 [ # # # # : 0 : nlrec%eival(:,:,:) = zero
# # ]
670 [ # # # # ]: 0 : nlrec%radii(:,:) = zero
671 [ # # # # ]: 0 : nlrec%pspinfo(:,:) = 0
672 : :
673 : : !--Get the exponential of the strength times the projectors overlap
674 : : ! of the non-local part of psp(hgh):
675 : : ! h_ij=strength; g_ij=ovelap => (exp(-h.g/temp/4p)-Identity).g^(-1)
676 : : ! And the diagonalisation of the projectors and associated eigenvectors
677 : 0 : call pspnl_hgh_rec(psps,tempe,nlrec,debug)
678 : :
679 [ # # ]: 0 : if(debug)then
680 [ # # ]: 0 : do jj=1,psps%npsp
681 : 0 : write(msg,'(a)')' Exponential matrices:'
682 : 0 : call wrtout(std_out,msg,'COLL')
683 [ # # ]: 0 : do ii=1,psps%mpsang
684 : : write(msg,'(a,i2,a,3f15.10,a,3f15.10,a,3f15.10)')&
685 : 0 : & 'angular moment',ii-1,ch10,&
686 : 0 : & nlrec%mat_exp_psp_nl(1,:,ii,jj),ch10,&
687 : 0 : & nlrec%mat_exp_psp_nl(2,:,ii,jj),ch10,&
688 : 0 : & nlrec%mat_exp_psp_nl(3,:,ii,jj)
689 : 0 : call wrtout(std_out,msg,'COLL')
690 : : end do
691 : : end do
692 : : end if
693 : :
694 : : !--Now it calculates the matrix of the exp(V_NL)
695 : 0 : call pspnl_operat_rec(nlrec,metrec,ngfftrec,debug)
696 : :
697 : : else !--Only local pseudo potentials
698 : 6 : nlrec%nlpsp = .False.
699 : 6 : nlrec%npsp = psps%npsp
700 [ - + - + ]: 6 : ABI_MALLOC(nlrec%mat_exp_psp_nl,(0,0,0,0))
701 [ - + - + ]: 6 : ABI_MALLOC(nlrec%pspinfo,(0,0))
702 [ - + - + ]: 6 : ABI_MALLOC(nlrec%radii,(0,0))
703 [ - + - + ]: 6 : ABI_MALLOC(nlrec%indlmn,(0,0,0))
704 [ - + - + ]: 6 : ABI_MALLOC(nlrec%projec,(0,0,0))
705 : : endif
706 : :
707 : 6 : end subroutine Init_nlpspRec
708 : : !!***
709 : :
710 : : !!****f* m_rec/CleanRec
711 : : !! NAME
712 : : !! CleanRec
713 : : !!
714 : : !! FUNCTION
715 : : !! Deallocate the pointers of rset<recursion_type>=Data type concerning recursion.
716 : : !!
717 : : !! INPUTS
718 : : !! rset<recursion_type>=Data type concerning recursion
719 : : !!
720 : : !! SIDE EFFECTS
721 : : !! All pointers are deallocated.
722 : : !!
723 : : !! SOURCE
724 : :
725 : 6 : subroutine CleanRec(rset)
726 : :
727 : : !Arguments ------------------------------------
728 : : ! scalars
729 : : type(recursion_type),intent(inout) :: rset
730 : : ! arrays
731 : : ! *********************************************************************
732 : :
733 : : ! @recursion_type
734 : :
735 [ + - ]: 6 : ABI_SFREE(rset%ZT_p)
736 [ + - ]: 6 : ABI_SFREE(rset%par%displs)
737 [ + - ]: 6 : ABI_SFREE(rset%par%vcount)
738 [ + - ]: 6 : ABI_SFREE(rset%nl%mat_exp_psp_nl)
739 [ - + ]: 6 : ABI_SFREE(rset%nl%eival)
740 [ - + ]: 6 : ABI_SFREE(rset%nl%eivec)
741 [ + - ]: 6 : ABI_SFREE(rset%nl%pspinfo)
742 [ + - ]: 6 : ABI_SFREE(rset%nl%radii)
743 [ + - ]: 6 : ABI_SFREE(rset%nl%indlmn)
744 [ + - ]: 6 : ABI_SFREE(rset%nl%projec)
745 [ + - ]: 6 : ABI_SFREE(rset%inf%gcart)
746 : :
747 : 6 : call pawfgr_destroy(rset%pawfgr)
748 : :
749 : : ! No is needed deallocate rset%mpi: it is a copy of mpi_enreg which
750 : : ! pointers are deallocated in gstate
751 : :
752 : : #ifdef HAVE_GPU_CUDA
753 : : call CleanRecGPU(rset%GPU,rset%load)
754 : : #endif
755 : :
756 : 6 : end subroutine CleanRec
757 : : !!***
758 : :
759 : : !!****f* m_rec/Calcnrec
760 : : !! NAME
761 : : !! Calcnrec
762 : : !!
763 : : !! FUNCTION
764 : : !! Calculate the new min_nrec.
765 : : !!
766 : : !! INPUTS
767 : : !! rset<recursion_type>=Data type concerning recursion
768 : : !! b2(:,:) recursion coefficients
769 : : !!
770 : : !! OUTPUT
771 : : !! rset%min_nrec is changed
772 : : !!
773 : : !! SOURCE
774 : :
775 : 12 : subroutine Calcnrec(rset,b2)
776 : :
777 : : !Arguments ------------------------------------
778 : : ! scalars
779 : : type(recursion_type),intent(inout) :: rset
780 : : ! arrays
781 : : real(dp), intent(in):: b2(0:rset%min_nrec,1:rset%par%ntranche)
782 : : !Local ----------------------------------------
783 : : ! scalars
784 : : integer :: kk,ii,jj,ierr,loc_nrec
785 : : character(len=500) :: msg
786 : : ! *********************************************************************
787 : : ! @recursion_type
788 : : ! @pawfgr_type
789 : 6 : kk = 1
790 : 6 : loc_nrec = rset%min_nrec
791 [ + + ]: 5190 : do ii=1,rset%par%ntranche
792 : : !--Use to lbound because b2 passed as argument
793 : : ! doesn't have the same bounds as in the calling
794 : : ! subroutine, the +1 because b2(lbound,ii)=1.
795 [ + - ]: 10368 : jj = lbound(b2,dim=1)+1
796 [ + + - + ]: 23631 : do while (b2(jj,ii)>tol10 .and. jj<=rset%min_nrec-1)
797 : 18441 : jj = jj+1
798 : 23625 : kk = max(jj,kk)
799 : : end do
800 : : enddo
801 : 6 : call xmpi_max(kk,rset%min_nrec,rset%mpi%comm_bandfft,ierr)
802 [ + - ]: 6 : rset%min_nrec = rset%min_nrec+1-lbound(b2,dim=1)
803 : :
804 : 6 : write(msg,'(a,i2.2,a9,i2.2)') ' -- nrec adjustement nrec=',loc_nrec,' => nrec=',rset%min_nrec
805 : 6 : call wrtout(std_out,msg,'COLL')
806 [ + - + + ]: 306 : write(msg,'(51a)')' ',('-',ii=1,50)
807 : 6 : call wrtout(std_out,msg,'COLL')
808 : :
809 : 6 : end subroutine Calcnrec
810 : : !!***
811 : :
812 : : !!****f* ABINIT/getngrec
813 : : !! NAME
814 : : !! getngrec
815 : : !!
816 : : !! FUNCTION
817 : : !! This routine computes the fft box for the recursion method, accordingly to the troncation radius.
818 : : !! It is quite similar to getng, but :
819 : : !! - there is no xboxcut and ecut consistency
820 : : !! - ngfft (the initial fft box) is the maximum fft box
821 : : !!
822 : : !! INPUTS
823 : : !! ngfft(18)=non truncated fft box
824 : : !! mgfft=maximum of ngfft(1:3)
825 : : !! inf_rmet=define the infinitesimal metric : rprimd*(transpose(rprimd)), divided by the number of discretisation point
826 : : !! recrcut=truncating
827 : : !! delta=to obtain radius of truncation
828 : : !!
829 : : !! OUTPUT
830 : : !! ngfftrec=truncated fft box
831 : : !! nfftrec= truncated nfft
832 : : !! tronc=True if truncation is made
833 : : !!
834 : : !! SIDE EFFECTS
835 : : !!
836 : : !! SOURCE
837 : :
838 : 6 : subroutine getngrec(ngfft,rmet,ngfftrec,nfftrec,recrcut,delta,tronc)
839 : :
840 : : !Arguments -------------------------------
841 : : !scalars
842 : : real(dp),intent(in) :: recrcut,delta
843 : : integer,intent(out) :: nfftrec
844 : : logical,intent(out) :: tronc
845 : : !arrays
846 : : integer,intent(in) :: ngfft(18)
847 : : real(dp),intent(in) :: rmet(3,3)
848 : : integer,intent(out) :: ngfftrec(18)
849 : :
850 : : !Local variables-------------------------------
851 : : !scalars
852 : : integer :: ii,iimin,index,jj,jjmin,kk,kkmin,largest_ngfftrec,maxpow11,maxpow2
853 : : integer :: maxpow3,maxpow5,maxpow7,mmsrch,plane
854 : : real(dp) :: dsm,dsp,dsqmin,rtroncat
855 : : !arrays
856 : : integer :: get_ngfftrec(3),maxsrch(3),minsrch(3)
857 : 6 : integer,allocatable :: iperm(:),srch(:)
858 : : real(dp) :: tsec(2)
859 : : real(dp) :: inf_rmet(3,3)
860 : : ! *************************************************************************
861 : :
862 : 6 : call timab(602,1,tsec)
863 : :
864 : 6 : ngfftrec(:) = ngfft(:)
865 : :
866 [ + + ]: 6 : if (recrcut>tol14) then !default value dtset%recrcut = zero means no troncation
867 : 1 : rtroncat = recrcut+delta
868 [ + + ]: 4 : get_ngfftrec(:)=1
869 : 1 : plane = 1
870 : :
871 [ + + ]: 4 : do ii=1,3
872 [ + + ]: 13 : inf_rmet(ii,:) = rmet(ii,:)/(real(ngfft(1:3)*ngfft(ii),dp))
873 : : end do
874 : :
875 : :
876 : : ! minimum value of ngfftrec
877 [ + + ]: 4 : do ii = 1,3
878 [ + - ]: 3 : ngfftrec(ii)=floor(2*rtroncat/sqrt(inf_rmet(ii,ii)))+1 !minimum value
879 [ - + ]: 4 : if(ngfftrec(ii)>=ngfft(ii))then
880 : 0 : ngfftrec(ii)=ngfft(ii)
881 : 0 : get_ngfftrec(ii)=0
882 : : end if
883 : : end do
884 : :
885 : :
886 [ + + + - ]: 4 : if(sum(get_ngfftrec)/=0)then
887 [ + + ]: 4 : largest_ngfftrec=maxval(ngfft(1:3))
888 : 1 : maxpow2=int(log(largest_ngfftrec+0.5d0)/log(two))
889 : 1 : maxpow3=int(log(largest_ngfftrec+0.5d0)/log(three))
890 : 1 : maxpow5=int(log(largest_ngfftrec+0.5d0)/log(five))
891 : 1 : maxpow7=0
892 : 1 : maxpow11=0
893 : 1 : mmsrch=(maxpow2+1)*(maxpow3+1)*(maxpow5+1)*(maxpow7+1)*(maxpow11+1)
894 [ - + + - : 3 : ABI_MALLOC(srch,(mmsrch))
+ - + - -
+ - + ]
895 [ - + + - : 2 : ABI_MALLOC(iperm,(mmsrch))
- + - + ]
896 : : ! Factors of 2
897 : 1 : srch(1)=1
898 [ + + ]: 4 : do ii=1,maxpow2
899 : 4 : srch(ii+1)=srch(ii)*2
900 : : end do
901 : : ! Factors of 3
902 : 1 : index=maxpow2+1
903 [ + - ]: 1 : if(maxpow3>0)then
904 [ + + ]: 3 : do ii=1,maxpow3
905 [ + - + + : 19 : srch(1+ii*index:(ii+1)*index)=3*srch(1+(ii-1)*index:ii*index)
+ + ]
906 : : end do
907 : : end if
908 : : ! Factors of 5
909 : 1 : index=(maxpow3+1)*index
910 [ + - ]: 1 : if(maxpow5>0)then
911 [ + + ]: 2 : do ii=1,maxpow5
912 [ + - + + : 26 : srch(1+ii*index:(ii+1)*index)=5*srch(1+(ii-1)*index:ii*index)
+ + ]
913 : : end do
914 : : end if
915 : : ! Factors of 7
916 : : index=(maxpow5+1)*index
917 : : if(maxpow7>0)then
918 : : do ii=1,maxpow7
919 : : srch(1+ii*index:(ii+1)*index)=7*srch(1+(ii-1)*index:ii*index)
920 : : end do
921 : : end if
922 : : ! Factors of 11
923 : 1 : index=(maxpow7+1)*index
924 : : if(maxpow11>0)then
925 : : do ii=1,maxpow11
926 : : srch(1+ii*index:(ii+1)*index)=11*srch(1+(ii-1)*index:ii*index)
927 : : end do
928 : : end if
929 : : ! srch is the set of allowed ngfftrec values
930 : :
931 : 1 : call sort_int(mmsrch,srch,iperm)
932 [ - + ]: 1 : ABI_FREE(iperm)
933 : :
934 [ + + ]: 4 : do ii=1,3
935 [ + - ]: 3 : if(get_ngfftrec(ii)==1)then
936 [ + - ]: 12 : do jj=1,mmsrch
937 [ + + ]: 12 : if(srch(jj)>=ngfftrec(ii))then
938 : 3 : minsrch(ii)=jj
939 : 3 : ngfftrec(ii)=srch(jj)
940 : 3 : exit
941 : : end if
942 : : end do
943 [ + - ]: 24 : do jj=minsrch(ii),mmsrch
944 [ + + ]: 24 : if(srch(jj)>ngfft(ii))then
945 : : ! since ngfftrec(ii)<ngfft(ii) for get_ngfftrec(ii)==1,
946 : : ! and srch(mmsrch)maxval(ngfft(1:3)),
947 : : ! that will appens in the range minsrch(ii),mmsrch
948 : 3 : maxsrch(ii)=jj-1
949 : 3 : exit
950 : : end if
951 : : end do
952 : : end if
953 : : ! since ngfft(ii) is in srch, we have here srch(maxsrch(ii))=ngfft(ii)
954 : : ! minsrch(ii), maxsrch(ii) is the range of index of srch in which we can
955 : : ! search ngfftrec(ii)
956 : :
957 [ - + ]: 4 : if(ngfftrec(ii)>=ngfft(ii))then
958 : 0 : ngfftrec(ii)=ngfft(ii)
959 : 0 : get_ngfftrec(ii)=0
960 : : end if
961 : : end do
962 : : end if
963 : :
964 : : ! verify that the entiere truncation sphere is in the fft box ;
965 : : ! but only in the dimension in which we do not consider the entiere fft box
966 [ + + + - ]: 4 : do while(sum(get_ngfftrec)/=0) !again...
967 : :
968 : : ! determining the minimum distance between 0 and the boundary
969 : : ! of the fft box
970 : : ! quite similar to the subroutine "bound", but only over the plane which
971 : : ! are not the whole fft box
972 : 1 : dsqmin=dsq_rec(ngfftrec(1)/2,-ngfftrec(2)/2,-ngfftrec(3)/2,inf_rmet)+0.01d0
973 : :
974 [ + - ]: 1 : if(get_ngfftrec(1)/=0)then
975 : : ! look at +/- g1 planes:
976 [ + + ]: 6 : do jj=-ngfftrec(2)/2,ngfftrec(2)/2
977 [ + + ]: 31 : do kk=-ngfftrec(3)/2,ngfftrec(3)/2
978 : 25 : dsp = dsq_rec(ngfftrec(1)/2, jj, kk,inf_rmet)
979 : 25 : dsm = dsq_rec( - ngfftrec(1)/2, jj, kk,inf_rmet)
980 [ + + ]: 25 : if (dsp<dsqmin) then
981 : 6 : dsqmin = dsp
982 : 6 : iimin = ngfftrec(1)/2
983 : 6 : jjmin = jj
984 : 6 : kkmin = kk
985 : 6 : plane=1
986 : : end if
987 [ - + ]: 30 : if (dsm<dsqmin) then
988 : 0 : dsqmin = dsm
989 : 0 : iimin = - ngfftrec(1)/2
990 : 0 : jjmin = jj
991 : 0 : kkmin = kk
992 : 0 : plane=1
993 : : end if
994 : : end do
995 : : end do
996 : : end if
997 : :
998 [ + - ]: 1 : if(get_ngfftrec(2)/=0)then
999 : : ! +/- g2 planes:
1000 [ + + ]: 6 : do ii=-ngfftrec(1)/2,ngfftrec(1)/2
1001 [ + + ]: 31 : do kk=-ngfftrec(3)/2,ngfftrec(3)/2
1002 : 25 : dsp = dsq_rec(ii,ngfftrec(2)/2,kk,inf_rmet)
1003 : 25 : dsm = dsq_rec(ii,-ngfftrec(2)/2,kk,inf_rmet)
1004 [ - + ]: 25 : if (dsp<dsqmin) then
1005 : 0 : dsqmin = dsp
1006 : 0 : iimin = ii
1007 : 0 : jjmin = ngfftrec(2)/2
1008 : 0 : kkmin = kk
1009 : 0 : plane=2
1010 : : end if
1011 [ - + ]: 30 : if (dsm<dsqmin) then
1012 : 0 : dsqmin = dsm
1013 : 0 : iimin = ii
1014 : 0 : jjmin = - ngfftrec(2)/2
1015 : 0 : kkmin = kk
1016 : 0 : plane=2
1017 : : end if
1018 : : end do
1019 : : end do
1020 : : end if
1021 : :
1022 [ + - ]: 1 : if(get_ngfftrec(3)/=0)then
1023 : : ! +/- g3 planes:
1024 [ + + ]: 6 : do ii=-ngfftrec(1)/2,ngfftrec(1)/2
1025 [ + + ]: 31 : do jj=-ngfftrec(2)/2,ngfftrec(2)/2
1026 : 25 : dsp = dsq_rec(ii,jj,ngfftrec(3)/2,inf_rmet)
1027 : 25 : dsm = dsq_rec(ii,jj,-ngfftrec(3)/2,inf_rmet)
1028 [ - + ]: 25 : if (dsp<dsqmin) then
1029 : 0 : dsqmin = dsp
1030 : 0 : iimin = ii
1031 : 0 : jjmin = jj
1032 : 0 : kkmin = ngfftrec(3)/2
1033 : 0 : plane=3
1034 : : end if
1035 [ - + ]: 30 : if (dsm<dsqmin) then
1036 : 0 : dsqmin = dsm
1037 : 0 : iimin = ii
1038 : 0 : jjmin = jj
1039 : 0 : kkmin = - ngfftrec(3)/2
1040 : 0 : plane=3
1041 : : end if
1042 : : end do
1043 : : end do
1044 : : end if
1045 : :
1046 [ - + ]: 1 : if(dsqmin>=rtroncat)then
1047 : : get_ngfftrec=0
1048 : : exit
1049 : : end if
1050 : :
1051 : : ! Fix nearest boundary
1052 [ # # ]: 0 : do ii=minsrch(plane),maxsrch(plane)
1053 [ # # ]: 0 : if (srch(ii)>=ngfftrec(plane)) then
1054 : : ! redefine ngfft(plane) to next higher choice
1055 : 0 : ngfftrec(plane)=srch(ii+1)
1056 : : ! verify if we cover the whole box
1057 [ # # ]: 0 : if(ngfftrec(plane)>=ngfft(plane))then
1058 : 0 : ngfftrec(plane)=ngfft(plane)
1059 : 0 : get_ngfftrec(plane)=0
1060 : : end if
1061 : : ! Exit the loop over ii
1062 : : exit
1063 : : end if
1064 : : end do
1065 : :
1066 : : end do
1067 : :
1068 [ + - ]: 1 : if (allocated(srch)) then
1069 : 1 : ABI_FREE(srch)
1070 : : end if
1071 : :
1072 : : ! if(mod(ngfftrec(1),16)/=0) then
1073 : : ! ngfftrec(1) = ngfftrec(1)+(16-mod(ngfftrec(1),16))
1074 : : ! ngfftrec(2:3) = ngfftrec(1)
1075 : : ! endif
1076 : :
1077 : 1 : ngfftrec(4)=2*(ngfftrec(1)/2)+1
1078 : 1 : ngfftrec(5)=2*(ngfftrec(2)/2)+1
1079 : 1 : ngfftrec(6)=ngfftrec(3)
1080 : :
1081 : : ! --algorithm
1082 : 1 : ngfftrec(7)=ngfft(7) ! to be improved for a better non-parallel algorithm - here it is automatically 401
1083 : 1 : ngfftrec(8)=ngfft(8)
1084 : :
1085 : : end if
1086 : :
1087 : : !--For now, recursion method doesn't use paralelism on FFT - which would require a great number of processors
1088 [ + + ]: 24 : nfftrec = product(ngfftrec(1:3))
1089 [ + + ]: 24 : ngfftrec(9:11) = (/0,1,0/) !--(/ paral, nproc, %me \)
1090 [ + + ]: 18 : ngfftrec(12:13) = ngfftrec(2:3) ! n2proc ! n3proc
1091 : :
1092 [ + + + + ]: 9 : tronc = all(ngfftrec(:3)/=ngfft(:3))
1093 : 6 : call timab(602,2,tsec)
1094 : :
1095 : : contains
1096 : :
1097 : 151 : function dsq_rec(ii,jj,kk,inf_rmet)
1098 : :
1099 : : real(dp) :: dsq_rec
1100 : : integer,intent(in) :: ii,jj,kk
1101 : : real(dp),intent(in) :: inf_rmet(3,3)
1102 : : dsq_rec=sqrt(&
1103 : : & inf_rmet(1,1)*dble(ii**2)&
1104 : : & +inf_rmet(2,2)*dble(jj**2)&
1105 : : & +inf_rmet(3,3)*dble(kk**2)&
1106 : : & +two*(inf_rmet(1,2)*dble(ii*jj)&
1107 : : & +inf_rmet(2,3)*dble(jj*kk)&
1108 : 151 : & +inf_rmet(3,1)*dble(kk*ii)))
1109 : 151 : end function dsq_rec
1110 : :
1111 : :
1112 : : end subroutine getngrec
1113 : : !!***
1114 : :
1115 : : !!****f* ABINIT/pspnl_operat_rec
1116 : : !! NAME
1117 : : !! pspnl_operat_rec
1118 : : !!
1119 : : !! FUNCTION
1120 : : !! It calculates the non-local projectors used in recursion for any psp non-local:
1121 : : !! The nl interaction in recursion is $$exp{-V_{NL}/beta}=\sum_A\sum_{lm}
1122 : : !! \sum{ij}Y_{lm}(\hat{r-R_A}')f^l_i(r-R_A)D^l_{i,j}Y_{lm}(\hat{r-R_A})f^l_j{r-R_A}$$
1123 : : !! where $D^_{i,j}$ is a matrix previously (see pspnl_operat_rec).
1124 : : !! In this routine the projectors $Y_{lm}(\hat{r-R_A}')f^l_i(r-R_A)$
1125 : : !! are calculated. So an array of dimensions
1126 : : !! rset%nl%projec(nfftrec,lmnmax,nlrec%npsp)
1127 : : !!
1128 : : !! INPUTS
1129 : : !! metrec<metricrec_type>=contains information concerning metric in
1130 : : !! recursion: grid_step, metric, infinitesimal volume
1131 : : !! ngfftrec(18)=is the ngfft grid (truncated if different from ngfft) of recursion
1132 : : !! debug=debug variable
1133 : : !!
1134 : : !!
1135 : : !! OUTPUT
1136 : : !! nlrec<nlpsprec_type>%projec= array containig the projectors on the real grid
1137 : : !! nlrec<nlpsprec_type>%intlen= integer linear size of the non-local grid
1138 : : !!
1139 : : !! SIDE EFFECTS
1140 : : !! nlrec<nlpsprec_type> data set of non-local pseudo for recursion
1141 : : !! The better Interaction length (Intlen) is also calculated and printed but
1142 : : !! recursion use intlen=ngfftrec/2
1143 : : !!
1144 : : !! NOTES
1145 : : !!
1146 : : !! SOURCE
1147 : :
1148 : :
1149 : 0 : subroutine pspnl_operat_rec(nlrec,metrec,ngfftrec,debug)
1150 : :
1151 : : !Arguments ------------------------------------
1152 : : !scalars
1153 : : logical,intent(in) :: debug
1154 : : type(metricrec_type),intent(in) ::metrec
1155 : : type(nlpsprec_type),intent(inout) :: nlrec
1156 : : !arrays
1157 : : integer,intent(in) :: ngfftrec(18)
1158 : : !Local variables-------------------------------
1159 : : integer :: ii,intlen
1160 : : integer :: iangol,ipsp,iproj
1161 : : integer :: mpsang,jj,kk,rr
1162 : : integer :: nfftrec
1163 : : integer :: ilmn,il,ilm,in,lmnmax
1164 : : real(dp) :: raggio,rloc,denom,step
1165 : : real(dp) :: delta_out,partial,err
1166 : : character(len=500) :: msg
1167 : : real(dp) :: part_sum(3)
1168 : : real(dp) :: ylmr_gr_dm(0,0,0)
1169 : 0 : real(dp),allocatable :: ylmr(:,:),proj_arr(:,:,:)
1170 : 0 : real(dp),allocatable :: radloc(:,:),nrm(:)
1171 : :
1172 : : ! *************************************************************************
1173 : :
1174 [ # # ]: 0 : if(debug)then
1175 [ # # # # ]: 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_operat_rec : enter '
1176 : 0 : call wrtout(std_out,msg,'PERS')
1177 : : end if
1178 : :
1179 : : !#####################################################################
1180 : : !--CALCULATE THE (SEMI-)MAXIMUM INTERVAL WHERE ALL THE PROJECTORS ARE
1181 : : !DIFFERENT TO ZERO.
1182 : : !--For any pseudo potential:
1183 : 0 : delta_out = zero
1184 : 0 : step = metrec%tr(1)*half !--Scanning step= grid step/2
1185 : :
1186 : :
1187 [ # # ]: 0 : do ipsp = 1, nlrec%npsp !--Loop on the pseudos
1188 : :
1189 : : ! --For any angular moment:
1190 [ # # # # ]: 0 : do iangol = 0,maxval(nlrec%indlmn(1,:,ipsp)) !--Loop on the angular moment
1191 : 0 : rloc = nlrec%radii(iangol+1,ipsp) !--Local radius
1192 : :
1193 : : ! --For any projector
1194 [ # # ]: 0 : do iproj = 1,nlrec%pspinfo(iangol+1,ipsp)
1195 : : ! --Starting point to searching when the projector goes to zero.
1196 : : ! this correspond to twice the radius wher the projector has its maximum
1197 : 0 : raggio = two*sqrt(real(-2+2*iproj+iangol,dp))*rloc
1198 : : ! --Caculate the gamma func at the denominator
1199 : 0 : call gamma_function(real(iangol+2*iproj,dp)-half,denom)
1200 : : ! --Find the zero
1201 : : ! --The following while cycle should be replaced by a bisection
1202 : : ! --method. Bucause this is calculated only 1 time it is not very
1203 : : ! important.
1204 : 0 : err = one
1205 : 0 : ii=0
1206 : : ! tolloc = 1.d0*abs(minval(nlrec%mat_exp_psp_nl(:nlrec%pspinfo(iangol+1,ipsp),:nlrec%pspinfo(iangol+1,ipsp),iangol+1,ipsp)))
1207 [ # # ]: 0 : do while(abs(err)>1.d-2)
1208 : 0 : raggio = raggio + step
1209 : 0 : err = project_prec(raggio,iproj,iangol,rloc)/sqrt(denom)
1210 : 0 : ii = ii+1
1211 : : end do
1212 : 0 : write(std_out,*)'local delta',raggio,ii
1213 [ # # # # : 0 : delta_out=maxval((/ delta_out,raggio /))
# # ]
1214 : : end do !end loop on projectors
1215 : :
1216 : : end do !enddo on angular moment
1217 : : end do !enddo on pseudos
1218 : :
1219 : : !--CALCULATE how many grid steps correspond to delta_out
1220 : 0 : intlen = int(delta_out/metrec%tr(1))
1221 : : !--I want that intlen is odd
1222 [ # # ]: 0 : if(mod(intlen,2)==0) intlen = intlen+1
1223 : :
1224 : 0 : write(msg,'(2a,i3,a)') ch10,' Interac. length of non-local psp(grid steps)=',intlen,ch10
1225 : 0 : call wrtout(std_out,msg,'COLL')
1226 : : !#####################################################################
1227 : :
1228 : : !--Initialisation
1229 [ # # ]: 0 : nfftrec = product(ngfftrec(1:3))
1230 : 0 : lmnmax = nlrec%lmnmax
1231 : 0 : intlen = ngfftrec(1)/2
1232 : 0 : nlrec%intlen = intlen !--Setted in recursion variables
1233 : :
1234 : : !#####################################################################
1235 : : !--CALCULATE E(q,q')
1236 : : !--Cration of the exponential*projectors*ylm matrix
1237 : :
1238 : : !--Initialisation
1239 [ # # # # : 0 : ABI_MALLOC(nlrec%projec,(nfftrec,lmnmax,nlrec%npsp))
# # # # #
# # # # #
# # # # #
# # # ]
1240 [ # # # # : 0 : nlrec%projec = zero
# # ]
1241 [ # # # # : 0 : ABI_MALLOC(radloc,(3,nfftrec))
# # # # #
# # # ]
1242 [ # # # # ]: 0 : radloc = zero
1243 [ # # # # : 0 : ABI_MALLOC(nrm,(nfftrec))
# # # # #
# # # ]
1244 [ # # ]: 0 : nrm = zero
1245 : :
1246 : : !--Loop on pseudo types
1247 [ # # ]: 0 : pseudodo: do ipsp = 1, nlrec%npsp
1248 : : ! --Control if the psp is non-local, else continue
1249 [ # # # # : 0 : if(all(nlrec%pspinfo(:,ipsp)==0)) cycle
# # ]
1250 : : ! --Vector which stores localy the upper part of symmetrical
1251 : : ! matrix of the exponential of the non-local operator
1252 [ # # ]: 0 : mpsang = maxval(nlrec%indlmn(1,:,ipsp))+1
1253 [ # # # # : 0 : ABI_MALLOC(proj_arr,(nfftrec,maxval(nlrec%pspinfo(:,ipsp)),mpsang))
# # # # #
# # # # #
# # # # #
# # # ]
1254 [ # # # # : 0 : ABI_MALLOC(ylmr,(mpsang*mpsang,nfftrec))
# # # # #
# # # # #
# # ]
1255 [ # # # # : 0 : proj_arr = zero
# # ]
1256 [ # # # # ]: 0 : ylmr = zero
1257 : :
1258 : : ! !debug
1259 : : ! write(std_out,*)'mpsang,proj num',mpsang,maxval(nlrec%pspinfo(:,ipsp))
1260 : : ! !enddebug
1261 : :
1262 : : ! --Calculate the projctors
1263 [ # # ]: 0 : do iangol = 0,mpsang-1
1264 : 0 : rloc = nlrec%radii(iangol+1,ipsp)
1265 [ # # ]: 0 : do iproj = 1,nlrec%pspinfo(iangol+1,ipsp)
1266 : 0 : call gamma_function(real(iangol+2*iproj,dp)-half,denom)
1267 : 0 : denom = one/sqrt(denom)
1268 [ # # ]: 0 : do ii = 0,ngfftrec(1)-1 !--3-loop on coordinates
1269 [ # # ]: 0 : do jj = 0,ngfftrec(2)-1
1270 [ # # ]: 0 : do kk = 0,ngfftrec(3)-1
1271 : : ! --Calculate the radii
1272 [ # # ]: 0 : part_sum(:) = real((/ ii,jj,kk /)-intlen,dp)*(metrec%tr)
1273 : 0 : rr = 1+ii+(jj+kk*ngfftrec(2))*ngfftrec(3)
1274 [ # # ]: 0 : radloc(:,rr) = part_sum
1275 [ # # ]: 0 : nrm(rr) = sqrt(sum(part_sum(:)**two))
1276 : 0 : partial = project_prec(nrm(rr),iproj,iangol,rloc)*denom
1277 [ # # ]: 0 : if(abs(partial)>tol12 ) proj_arr(rr,iproj,iangol+1) = partial
1278 : : end do
1279 : : end do
1280 : : end do !--End 3-loop on coordinates
1281 : : end do
1282 : : end do
1283 : :
1284 : :
1285 : : ! -------------------------------------------------------------
1286 : : ! --Calculate the spherical harmonics (Verified: it works well)
1287 : 0 : call initylmr(mpsang,1,nfftrec,nrm(:),1,radloc(:,:),ylmr(:,:),ylmr_gr_dm)
1288 : : ! -------------------------------------------------------------
1289 : :
1290 : :
1291 [ # # ]: 0 : do ilmn = 1,lmnmax
1292 : 0 : ilm = nlrec%indlmn(4,ilmn,ipsp)
1293 : 0 : il = nlrec%indlmn(1,ilmn,ipsp)+1
1294 : 0 : in = nlrec%indlmn(3,ilmn,ipsp)
1295 : 0 : write(msg,'(2a,i3,2i2)')ch10,'lm,l,n',ilm,il,in
1296 : 0 : call wrtout(std_out,msg,'COLL')
1297 : :
1298 [ # # ]: 0 : nlrec%projec(:,ilmn,ipsp) = ylmr(ilm,:)*proj_arr(:,in,il)
1299 : : end do
1300 : :
1301 [ # # ]: 0 : ABI_FREE(ylmr)
1302 [ # # ]: 0 : ABI_FREE(proj_arr)
1303 : : end do pseudodo !--end loop on pseudo types
1304 : :
1305 : :
1306 [ # # ]: 0 : ABI_FREE(radloc)
1307 [ # # ]: 0 : ABI_FREE(nrm)
1308 : :
1309 [ # # ]: 0 : if(debug)then
1310 [ # # # # ]: 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_operat_rec : exit '
1311 : 0 : call wrtout(std_out,msg,'PERS')
1312 : : end if
1313 : :
1314 : : contains
1315 : :
1316 : 0 : function project_prec(raggio,iproj,iangol,rloc)
1317 : : !--Analytical expression of the projectors in hgh-pspeudopotential
1318 : : !--The gamma function at denominator is missing
1319 : : real(dp) :: project_prec
1320 : : integer,intent(in) :: iproj,iangol
1321 : : real(dp),intent(in) :: raggio,rloc
1322 : :
1323 : : project_prec=sqrt2*(raggio/rloc)**real((iangol+2*(iproj-1)),dp)*&
1324 : 0 : & exp(-((raggio/rloc)**two)*half)/rloc**onehalf
1325 : 0 : end function project_prec
1326 : :
1327 : : end subroutine pspnl_operat_rec
1328 : : !!***
1329 : :
1330 : : !!****f* ABINIT/pspnl_hgh_rec
1331 : : !! NAME
1332 : : !! pspnl_hgh_rec
1333 : : !!
1334 : : !! FUNCTION
1335 : : !! This routine computes the matrices S_kk'^{l,A} of the projectors
1336 : : !! (it is the exponential of the overlap matrix). It coorresponds to the matrix:
1337 : : !! $$\left[(U_l)^{-1}*Exp(-temperature D_l )*U_l* (g_l)^{-1} -Identity\right]_kk'
1338 : : !! where (U_l)^-1* D_l* U_l = h^lg_l.
1339 : : !! $g_l = <f^l_k|f^l_{k'}>$ is the overlap matrix between projectors
1340 : : !! and $h^l_{kk'}$ is the strength matrix of the projectors.
1341 : : !! It calulates also the strength eigenvalues and eigenvectors of $h$,
1342 : : !! used in the calculus of non-local energy
1343 : : !!
1344 : : !! INPUTS
1345 : : !! temperature=4*rtrotter/beta=4*rtrotter*tsmear: the effective temp. in recursion
1346 : : !! psps <type(pseudopotential_type)>=variables related to pseudo-potentials
1347 : : !! debug=debug variable
1348 : : !!
1349 : : !! OUTPUT
1350 : : !!
1351 : : !! nlrec%mat_exp_psp_nl=the matrix of the exponential of the projectors:
1352 : : !! for any psp, for any angular moment:
1353 : : !! h_ij=strength; g_ij=ovelap => exp(-h.g/temp/4p).g^(-1)
1354 : : !! nlrec%radii=Local radii of nl psp
1355 : : !! nlrec%pspinfo(:,:) for any typat: (momang,typat)=number of projectors
1356 : : !! nlrec%eival(:,:,:) for any psp, any momang, the eigenvalues of the
1357 : : !! strength matrix H: eival(:,mang,psp)
1358 : : !! nlrec%eivec(:,:,:,:)for any psp, any momang, the eigenvectors of the
1359 : : !! strength matrix H: eivec(:,:,mang,psp)
1360 : : !!
1361 : : !! SIDE EFFECTS
1362 : : !!
1363 : : !! SOURCE
1364 : :
1365 : 0 : subroutine pspnl_hgh_rec(psps,temperature,nlrec,debug)
1366 : :
1367 : : use m_linalg_interfaces
1368 : :
1369 : : !Arguments -----------------------------------
1370 : : !scalars
1371 : : real(dp),intent(in) :: temperature
1372 : : logical,intent(in) :: debug
1373 : : type(pseudopotential_type),intent(in) :: psps
1374 : : type(nlpsprec_type),intent(inout) :: nlrec
1375 : : !arrays
1376 : : !Local variables-------------------------------
1377 : : !scalars
1378 : : integer,parameter :: maxsize=3
1379 : : integer,parameter :: lwork=(1+32)*maxsize
1380 : : integer :: iangol,ipseudo,info,nproj
1381 : : integer :: g_mat_size,ii,nproj2
1382 : : real(dp) :: denom_1,denom_2,tot_proj
1383 : : character(len=500) :: msg
1384 : : !arrays
1385 : : integer :: ipvt(1)
1386 : : !real(dp) :: rwork(2*maxsize)
1387 : : real(dp) :: h_mat_init(3,3), rework(lwork)
1388 : 0 : real(dp), allocatable :: g_mat(:,:),h_mat(:,:),eig_val_h(:)
1389 : 0 : real(dp), allocatable :: identity(:,:),inv_g_mat(:,:),u_mat(:,:)
1390 : 0 : complex(dp),allocatable :: hg_mat(:,:)
1391 : : ! *************************************************************************
1392 : :
1393 [ # # ]: 0 : if(debug)then
1394 [ # # # # ]: 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_hgh_rec : enter '
1395 : 0 : call wrtout(std_out,msg,'COLL')
1396 : : end if
1397 : :
1398 : : !--For any pseudo potential:
1399 [ # # ]: 0 : do ipseudo = 1, psps%npsp !--Loop on the pseudos
1400 [ # # # # ]: 0 : write(msg,'(a,80a)')' pseudo file',('-',ii=1,10)
1401 : 0 : call wrtout(std_out,msg,'COLL')
1402 : 0 : write(msg,'(a)') psps%filpsp(ipseudo)
1403 : 0 : call wrtout(std_out,msg,'COLL')
1404 : :
1405 : : ! --For any angular moment:
1406 [ # # ]: 0 : do iangol = 0,psps%mpsang-1 !--Loop on the angular moment
1407 : :
1408 : : ! --Local radius
1409 : 0 : nlrec%radii(iangol+1,ipseudo) = psps%gth_params%psppar(iangol+1,0,ipseudo)
1410 : :
1411 : : ! --Strenghts of non-local projectors (matrix h)
1412 : : ! --Diagonal part:
1413 : 0 : h_mat_init = zero
1414 : 0 : h_mat_init(1,1) = psps%gth_params%psppar(iangol+1,1,ipseudo)
1415 : 0 : h_mat_init(2,2) = psps%gth_params%psppar(iangol+1,2,ipseudo)
1416 : 0 : h_mat_init(3,3) = psps%gth_params%psppar(iangol+1,3,ipseudo)
1417 : : ! --Out-diagonal part
1418 : : ! --Depending on angular moment the projectors
1419 : : ! strength is calculated differently
1420 : 0 : select case(iangol)
1421 : : case(0)
1422 : 0 : h_mat_init(1,2) = -half*sqrt(3.d0/5.d0)*h_mat_init(2,2)
1423 : 0 : h_mat_init(1,3) = half*sqrt(5.d0/21.d0)*h_mat_init(3,3)
1424 : 0 : h_mat_init(2,3) = -half*sqrt(100.d0/63.d0)*h_mat_init(3,3)
1425 : : case(1)
1426 : 0 : h_mat_init(1,2) = -half*sqrt(5.d0/7.d0)*h_mat_init(2,2)
1427 : 0 : h_mat_init(1,3) = sixth*sqrt(35.d0/11.d0)*h_mat_init(3,3)
1428 : 0 : h_mat_init(2,3) = -14.d0/six/sqrt(11.d0) *h_mat_init(3,3)
1429 : : case(2)
1430 : 0 : h_mat_init(1,2) = -half*sqrt(7.d0/9.d0)*h_mat_init(2,2)
1431 : 0 : h_mat_init(1,3) = half*sqrt(63.d0/143.d0)*h_mat_init(3,3)
1432 : 0 : h_mat_init(2,3) = -nine/sqrt(143.d0)*h_mat_init(3,3)
1433 : : case(3)
1434 : 0 : h_mat_init(1,2) = zero; h_mat_init(1,3) = zero; h_mat_init(2,3) = zero;
1435 : : case default
1436 : 0 : write(msg,'(a)')' error angular: moment component'
1437 [ # # # # : 0 : call wrtout(std_out,msg,'COLL')
# ]
1438 : : end select
1439 : :
1440 : :
1441 : :
1442 : : ! --Real dimensions of projectors.
1443 [ # # # # ]: 0 : g_mat_size = count(abs((/ (h_mat_init(ii,ii),ii=1,3) /))>1.d-8)
1444 : 0 : nlrec%pspinfo(iangol+1,ipseudo) = g_mat_size
1445 : 0 : write(msg,'(a,i2,a,i2)')' ang. moment=',iangol,', N projectors=',g_mat_size
1446 : 0 : call wrtout(std_out,msg,'COLL')
1447 [ # # ]: 0 : if (g_mat_size>0) then
1448 : : ! --Identity matrix
1449 [ # # # # : 0 : ABI_MALLOC(identity,(g_mat_size,g_mat_size))
# # # # #
# # # # #
# # ]
1450 : 0 : call set2unit(identity)
1451 : : ! identity = zero
1452 : : ! identity(:,1) = one
1453 : : ! identity(:,:) = cshift(array=identity,shift=(/ (-ii,ii=0,g_mat_size) /), dim=2 )
1454 : :
1455 : :
1456 : : ! ############## CALCULOUS OF THE EIGEN_SPACE OF THE PROJECTORS STRENGTHS ##################
1457 : : ! --Inverse of the matrix h
1458 [ # # # # : 0 : ABI_MALLOC(eig_val_h,(g_mat_size))
# # # # #
# # # ]
1459 [ # # # # : 0 : ABI_MALLOC(u_mat,(g_mat_size,g_mat_size))
# # # # #
# # # ]
1460 : : ! --u-mat will contain the eigenvectors of h_mat_init
1461 [ # # # # : 0 : u_mat = h_mat_init(:g_mat_size,:g_mat_size)
# # ]
1462 : :
1463 : : ! write(std_out,*)'hmat_init'
1464 : : ! do ii=1,g_mat_size
1465 : : ! write(std_out,*)h_mat_init(ii,:)
1466 : : ! end do
1467 : 0 : call DSYEV('v','u',g_mat_size,u_mat,g_mat_size,eig_val_h,rework,lwork,info)
1468 : :
1469 : : ! --THE DIAGONAL MATRIX IS GIVEN BY D=U^t.H.U
1470 : : ! (eival=transpose(eivec).h_mat_init.eivec)
1471 : 0 : write(msg,'(a,3d10.3)')' eigenvalues=',eig_val_h
1472 : 0 : call wrtout(std_out,msg,'COLL')
1473 : : ! write(std_out,*)'autovec';write(std_out,*)u_mat
1474 : :
1475 [ # # ]: 0 : nlrec%eival(:g_mat_size,1+iangol,ipseudo) = eig_val_h
1476 [ # # # # ]: 0 : nlrec%eivec(:g_mat_size,:g_mat_size,1+iangol,ipseudo) = u_mat
1477 [ # # ]: 0 : ABI_FREE(eig_val_h)
1478 [ # # ]: 0 : ABI_FREE(u_mat)
1479 : :
1480 : : ! ##########END CALCULOUS OF THE EIGEN_SPACE OF THE PROJECTORS STRENGTHS ##################
1481 : :
1482 [ # # # # : 0 : ABI_MALLOC(g_mat,(g_mat_size,g_mat_size))
# # # # #
# # # # #
# # # # ]
1483 [ # # # # : 0 : ABI_MALLOC(inv_g_mat,(g_mat_size,g_mat_size))
# # # # #
# # # #
# ]
1484 [ # # # # : 0 : ABI_MALLOC(h_mat,(g_mat_size,g_mat_size))
# # # # #
# # # #
# ]
1485 [ # # # # : 0 : ABI_MALLOC(hg_mat,(g_mat_size,g_mat_size))
# # # # #
# # # # #
# # # # ]
1486 : :
1487 [ # # # # ]: 0 : g_mat(:,:) = one
1488 [ # # # # ]: 0 : h_mat(:,:) = zero
1489 [ # # # # ]: 0 : h_mat(:,:) = h_mat_init(:g_mat_size,:g_mat_size)
1490 : :
1491 : : ! -------------------------------------------------------
1492 : : ! --Matrix of the overlap between projetors (matrix g)
1493 : : ! and the h matrix of strength
1494 [ # # ]: 0 : do nproj = 1,g_mat_size-1
1495 [ # # ]: 0 : do nproj2 = 1+nproj,g_mat_size
1496 : : tot_proj = zero
1497 : : ! --Analytic value of overlap
1498 : : ! g_ij=Gamma[-1/2+i+j+l]/Sqrt(Gamma[-1/2+i+iangol]*Gamma[-1/2+j+iangol])
1499 : 0 : call gamma_function(-half+real(nproj+nproj2+iangol,dp),tot_proj)
1500 : 0 : call gamma_function(-half+real(iangol+2*nproj,dp),denom_1)
1501 : 0 : call gamma_function(-half+real(iangol+2*nproj2,dp),denom_2)
1502 : :
1503 : 0 : g_mat(nproj,nproj2) = tot_proj/sqrt(denom_1*denom_2)
1504 : 0 : g_mat(nproj2,nproj) = g_mat(nproj,nproj2)
1505 : :
1506 : 0 : h_mat(nproj,nproj2) = h_mat_init(nproj,nproj2)
1507 : 0 : h_mat(nproj2,nproj) = h_mat_init(nproj,nproj2)
1508 : : end do
1509 : : end do
1510 : :
1511 : : ! --Inverse of the overlap matrix g
1512 [ # # # # : 0 : inv_g_mat = g_mat
# # ]
1513 : 0 : call DGETRF(g_mat_size,g_mat_size,inv_g_mat,g_mat_size,ipvt,info)
1514 : 0 : call DGETRI(g_mat_size,inv_g_mat,g_mat_size,ipvt,rework,lwork,info)
1515 : :
1516 : :
1517 : : ! -----------------------------------------------------------
1518 : : ! --Now it calculates the exponential of the matrix h.g
1519 [ # # # # : 0 : hg_mat = matmul(h_mat,g_mat)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
1520 : :
1521 : 0 : call exp_mat(hg_mat,g_mat_size,-one/temperature)
1522 : :
1523 : : ! --(exp(h.g)-Identity).(g^-1)
1524 [ # # # # : 0 : hg_mat = hg_mat-identity(:,:)
# # ]
1525 : :
1526 : :
1527 : : ! --results on output
1528 [ # # # # : 0 : nlrec%mat_exp_psp_nl(:g_mat_size,:g_mat_size,1+iangol,ipseudo) = matmul(real(hg_mat,dp),inv_g_mat)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
1529 : :
1530 : : ! write(std_out,*) nlrec%mat_exp_psp_nl(:g_mat_size,:g_mat_size,1+iangol,ipseudo)
1531 : :
1532 [ # # ]: 0 : ABI_FREE(g_mat)
1533 [ # # ]: 0 : ABI_FREE(hg_mat)
1534 [ # # ]: 0 : ABI_FREE(h_mat)
1535 [ # # ]: 0 : ABI_FREE(inv_g_mat)
1536 [ # # ]: 0 : ABI_FREE(identity)
1537 : : end if
1538 : :
1539 : : end do !enddo on angular moment
1540 : : end do !enddo on pseudos
1541 : :
1542 : :
1543 [ # # ]: 0 : if(debug)then
1544 [ # # # # ]: 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_hgh_rec : exit '
1545 : 0 : call wrtout(std_out,msg,'COLL')
1546 : : end if
1547 : :
1548 [ # # # # : 0 : end subroutine pspnl_hgh_rec
# # # # #
# # # ]
1549 : : !!***
1550 : :
1551 [ # # ]: 0 : end module m_rec
1552 : : !!***
|
