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