Line data Source code
1 : !!****m* ABINIT/m_vtorhorec
2 : !! NAME
3 : !! m_vtorhorec
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (SLeroux, MMancini).
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_vtorhorec
22 :
23 : use defs_basis
24 : use defs_rectypes
25 : use m_xmpi
26 : use m_pretty_rec
27 : use m_errors
28 : use m_abicore
29 : use m_per_cond
30 : use m_dtset
31 :
32 : use defs_datatypes, only : pseudopotential_type
33 : use defs_abitypes, only : MPI_type
34 : use m_time, only : timein, timab
35 : use m_rec, only : Calcnrec, init_nlpsprec, cpu_distribution
36 : use m_rec_tools, only : reshape_pot, trottersum, get_pt0_pt1
37 : use m_spacepar, only : symrhg
38 : use m_fourier_interpol, only : transgrid
39 : use m_fft, only : fourdp
40 :
41 : #ifdef HAVE_GPU_CUDA
42 : use m_gpu_toolbox
43 : use m_hidecudarec
44 : use m_xredistribute
45 : #endif
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : public :: vtorhorec
53 : public :: first_rec
54 : !!***
55 :
56 : contains
57 : !!***
58 :
59 : !!****f* ABINIT/vtorhorec
60 : !! NAME
61 : !! vtorhorec
62 : !!
63 : !! FUNCTION
64 : !! This routine computes the new density from a fixed potential (vtrial)
65 : !! using a recursion method
66 : !!
67 : !! INPUTS
68 : !! deltastep= if 0 the iteration step is equal to dtset%nstep
69 : !! initialized= if 0 the initialization of the gstate run is not yet finished
70 : !! operator (ground-state symmetries)
71 : !! dtset <type(dataset_type)>=all input variables for this dataset
72 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
73 : !! nfftf=(effective) number of FFT grid points (for this processor)
74 : !! nspden=number of spin-density components
75 : !! nsppol=1 for unpolarized, 2 for spin-polarized
76 : !! nsym=number of symmetry elements in space group
77 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
78 : !! vtrial(nfft,nspden)=INPUT Vtrial(r).
79 : !! rset <type(recursion_type)> all variables for recursion
80 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
81 : !! gprimd(3,3)=dimensional primitive translations in reciprocal space
82 : !!
83 : !! OUTPUT
84 : !! ek=kinetic energy part of total energy.
85 : !! enlx=nonlocal psp + potential Fock ACE part of total energy.
86 : !! entropy=entropy due to the occupation number smearing (if metal)
87 : !! e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree)
88 : !! fermie=fermi energy (Hartree)
89 : !! grnl(3*natom)=stores grads of nonlocal energy wrt length scales
90 : !! (3x3 tensor) and grads wrt atomic coordinates (3*natom)
91 : !!
92 : !! SIDE EFFECTS
93 : !! rhog(2,nfft)=array for Fourier transform of electron density
94 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
95 : !! rset%efermi= fermi energy
96 : !!
97 : !! NOTES
98 : !! at this time :
99 : !! - grnl in not implemented
100 : !! - symetrie usage not implemented (irrzon not used and nsym should be 1)
101 : !! - spin-polarized not implemented (nsppol must be 1, nspden ?)
102 : !! - phnons used only in symrhg
103 : !! - need a rectangular box (ngfft(1)=ngfft(2)=ngfft(3))
104 : !!
105 : !! SOURCE
106 :
107 38 : subroutine vtorhorec(dtset,&
108 : & ek,enlx,entropy,e_eigenvalues,fermie,&
109 76 : & grnl,initialized,irrzon,nfftf,phnons,&
110 76 : & rhog, rhor, vtrial,rset,deltastep,rprimd,gprimd)
111 :
112 : !Arguments -------------------------------
113 : !scalars
114 : integer,intent(in) :: initialized
115 : integer,intent(in) :: nfftf,deltastep
116 : real(dp),intent(out) :: e_eigenvalues,ek,enlx,entropy,fermie
117 : type(dataset_type),intent(in) :: dtset
118 : type(recursion_type),intent(inout) :: rset
119 : !arrays
120 : integer, intent(in) :: irrzon(:,:,:)
121 : real(dp),intent(in) :: rprimd(3,3),gprimd(3,3)
122 : real(dp),intent(in) :: phnons(:,:,:)
123 : real(dp),intent(in) :: vtrial(:,:)
124 : real(dp),intent(inout) :: rhog(:,:)
125 : real(dp),intent(out) :: grnl(:) !vz_i
126 : real(dp),intent(inout) :: rhor(:,:) !vz_i
127 :
128 : !Local variables-------------------------------
129 : !scalars
130 : integer :: nfftrec, dim_entro,ii1,jj1,kk1
131 : integer :: ierr,ii,ilmn,ipsp,dim_trott
132 : integer :: ipoint,ipointlocal,jj,kk,irec
133 : integer :: n1,n2,n3,n4,n_pt_integ_entropy
134 : integer :: nrec,iatom,min_pt,max_pt,swt_tm
135 : integer :: get_K_S_G
136 : integer :: tim_fourdp,trotter
137 : real(dp),parameter :: perc_vmin=one
138 : real(dp) :: beta,drho,drhomax
139 : real(dp) :: entropy1,entropy2,entropy3,entropy4
140 : real(dp) :: entropylocal,entropylocal1,entropylocal2
141 : real(dp) :: entropylocal3,entropylocal4,gran_pot,gran_pot1,gran_pot2
142 : real(dp) :: gran_pot3,gran_pot4,gran_pot_local,gran_pot_local1
143 : real(dp) :: gran_pot_local2,gran_pot_local3,gran_pot_local4
144 : real(dp) :: inf_ucvol,intrhov,factor
145 : real(dp) :: nelect,potmin,rtrotter,toldrho,tolrec,tsmear
146 : real(dp) :: xmax,nlpotmin,ratio1,ratio2,ratio4,ratio8
147 38 : type(recparall_type) :: recpar
148 : character(len=500) :: msg
149 : !arrays
150 : integer :: ngfftrec(18),trasl(3)
151 38 : integer,pointer :: gcart_loc(:,:)
152 38 : integer,allocatable :: bufsize(:), bufdispl(:)
153 : real(dp) :: tsec(2),tsec2(2)
154 76 : real(dp) :: exppot(0:dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)-1)
155 76 : real(dp),target :: rholocal(1:rset%par%ntranche)
156 76 : real(dp),target :: alocal(0:rset%min_nrec,1:rset%par%ntranche)
157 76 : real(dp),target :: b2local(0:rset%min_nrec,1:rset%par%ntranche)
158 38 : real(dp),pointer :: rho_wrk(:)
159 38 : real(dp),pointer :: a_wrk(:,:),b2_wrk(:,:)
160 38 : real(dp),allocatable :: rholocal_f(:), rholoc_2(:)
161 38 : real(dp),allocatable :: rhogf(:,:),rhogc(:,:),aloc_copy(:,:),b2loc_copy(:,:)
162 38 : real(dp),allocatable :: gran_pot_v_f(:,:),gran_pot_v_c(:,:),gran_pot_v_2(:,:)
163 38 : real(dp),allocatable :: entropy_v_f(:,:),entropy_v_c(:,:),entropy_v_2(:,:)
164 38 : real(dp),allocatable :: ablocal_1(:,:,:),ablocal_2(:,:,:),ablocal_f(:,:,:)
165 38 : real(dp),allocatable :: exppotloc(:)
166 38 : real(dp),allocatable :: projec(:,:,:,:,:)
167 : #if defined HAVE_GPU_CUDA
168 : integer :: max_rec
169 : integer,allocatable :: vcount_0(:), displs_0(:)
170 : integer,allocatable :: vcount_1(:), displs_1(:)
171 : real(dp),allocatable,target :: rho_hyb(:)
172 : real(dp),allocatable,target :: a_hyb(:,:),b2_hyb(:,:)
173 : real(cudap),allocatable :: an_dev(:,:)
174 : real(cudap),allocatable :: bn2_dev(:,:)
175 : #endif
176 :
177 : ! *********************************************************************
178 38 : if(rset%debug)then
179 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' vtorhorec : enter '
180 0 : call wrtout(std_out,msg,'PERS')
181 : end if
182 :
183 38 : call timab(21,1,tsec)
184 38 : call timab(600,1,tsec2)
185 : !##################################################################################################
186 : !!--Initalization in the FIRST time in VTORHOREC is made in SCFCV by FIRST_REC routine
187 : !!--Parameters for the recursion method AND Initialisation
188 :
189 38 : trotter = dtset%recptrott !--trotter parameter
190 38 : nelect = dtset%nelect !--number of electrons
191 :
192 38 : toldrho = dtset%rectolden !--tollerance for density
193 38 : tolrec = toldrho*1.d-2 !--tollerance for local density
194 :
195 38 : tsmear = dtset%tsmear !--temperature
196 38 : beta = one/tsmear !--inverse of temperature
197 :
198 38 : factor = real(dtset%recgratio**3*100*rset%mpi%nproc,dp)/real(nfftf,dp)
199 :
200 : !--Assignation of the rset variable:
201 38 : nrec = rset%min_nrec
202 38 : nfftrec = rset%nfftrec
203 722 : ngfftrec = rset%ngfftrec
204 38 : inf_ucvol = rset%inf%ucvol
205 :
206 38 : min_pt = rset%par%displs(rset%mpi%me)+1
207 38 : max_pt = min_pt+rset%par%vcount(rset%mpi%me)-1
208 :
209 : !--In the last self-constistent loop or if density is converged the
210 : !thermodynamics quantities are calculated
211 38 : get_K_S_G = 0; if(deltastep==0 .or. rset%quitrec/=0) get_K_S_G = 1;
212 :
213 : !--Rewriting the trotter parameter
214 38 : rtrotter = max(half,real(trotter,dp))
215 38 : dim_trott = max(0,2*trotter-1)
216 :
217 : !--Variables Optimisation
218 38 : ratio1 = beta/rtrotter
219 38 : ratio2 = ratio1/two
220 38 : ratio4 = ratio1/four
221 38 : ratio8 = ratio1/eight
222 :
223 : !--Integration points for entropy
224 38 : n_pt_integ_entropy = max(25,dtset%recnpath)
225 :
226 : !-- energies non-local: at day not implemented
227 38 : enlx = zero
228 467 : grnl = zero
229 : !jmb
230 38 : ek = zero
231 :
232 : !--only a copy of ngfft(1:3) and nfft (to purge!!)
233 38 : n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3)
234 38 : n4 = n3
235 :
236 : !--time switch to measure gpu-cpu syncrhonisation
237 38 : swt_tm = 0 ; !no gpu initally
238 :
239 65702 : exppot = zero
240 38 : nullify(gcart_loc)
241 :
242 38 : if(dtset%rectesteg==1)then
243 : ! --Free electron gas case
244 5187 : exppot = one
245 : else
246 35 : if(.not.(rset%nl%nlpsp)) then
247 : ! --Local case
248 : ! --COMPUTATION OF exp( -beta*pot/(4*rtrotter))
249 35 : ABI_MALLOC(gcart_loc,(0,0))
250 : gcart_loc = 0
251 60515 : exppot = exp( -(ratio4*vtrial(:,1)))
252 : else
253 : ! --Non-Local case
254 : ! --COMPUTATION OF exp(-beta*pot/(8*rtrotter))
255 0 : exppot = exp( -(ratio8*vtrial(:,1)))
256 :
257 0 : ABI_MALLOC(gcart_loc,(3,dtset%natom))
258 0 : gcart_loc = rset%inf%gcart
259 0 : ABI_MALLOC(projec,(0:ngfftrec(1)-1,0:ngfftrec(2)-1,0:ngfftrec(3)-1,rset%nl%lmnmax,dtset%natom))
260 0 : projec = zero
261 :
262 0 : if(.not.(rset%tronc)) then
263 0 : do iatom =1, dtset%natom
264 0 : ipsp = dtset%typat(iatom)
265 0 : do ilmn = 1,rset%nl%lmnmax
266 0 : projec(:,:,:,ilmn,iatom) = reshape(rset%nl%projec(:,ilmn,ipsp),shape=shape(projec(:,:,:,1,1)))
267 0 : do ii=1,3
268 0 : projec(:,:,:,ilmn,iatom) = cshift(projec(:,:,:,ilmn,iatom),shift=ngfftrec(ii)/2-gcart_loc(ii,iatom),dim=ii)
269 : end do
270 : end do
271 : end do
272 : end if
273 : end if
274 : end if
275 :
276 : !###################################################################################
277 : !MAIN LOOP
278 :
279 955622 : rholocal = zero; alocal = zero; b2local = zero
280 38 : ipointlocal = 1
281 :
282 : !--Allocation: if hybrid calculation is done then I have to use
283 : !balanced work on devices.
284 38 : nullify(rho_wrk,a_wrk,b2_wrk)
285 :
286 38 : if(rset%load == 1)then
287 : #ifdef HAVE_GPU_CUDA
288 : ABI_MALLOC(rho_hyb,(1:rset%GPU%par%npt))
289 : ABI_MALLOC(a_hyb,(0:nrec,1:rset%GPU%par%npt))
290 : ABI_MALLOC(b2_hyb,(0:nrec,1:rset%GPU%par%npt))
291 : rho_hyb = zero; a_hyb = zero; b2_hyb = zero
292 : rho_wrk => rho_hyb
293 : a_wrk => a_hyb
294 : b2_wrk => b2_hyb
295 : recpar = rset%GPU%par
296 : #endif
297 : else
298 38 : rho_wrk => rholocal
299 38 : a_wrk => alocal
300 38 : b2_wrk => b2local
301 38 : recpar = rset%par
302 : end if
303 :
304 : !#if defined HAVE_GPU_CUDA
305 : !if(rset%debug)then
306 : !write (std_out,*) 'rset%recGPU%nptrec ',rset%recGPU%nptrec
307 : !write (std_out,*) 'rset%gpudevice ',rset%gpudevice
308 : !write (std_out,*) 'rset%ngfftrec ',rset%ngfftrec(1:3)
309 : !write (std_out,*) 'rset%min_nrec ',rset%min_nrec
310 : !write (std_out,*) 'rset%par%ntranche ',rset%par%ntranche
311 : !write (std_out,*) 'rset%par%min_pt ',rset%par%min_pt,min_pt
312 : !write (std_out,*) 'rset%par%max_pt ',rset%par%max_pt,max_pt
313 : !write (std_out,*) 'pt0 ',rset%par%pt0%x,rset%par%pt0%y,rset%par%pt0%z
314 : !write (std_out,*) 'pt1 ',rset%par%pt1%x,rset%par%pt1%y,rset%par%pt1%z
315 : !end if
316 : !#endif
317 :
318 38 : if(rset%gpudevice>=0) then
319 : #if defined HAVE_GPU_CUDA
320 : swt_tm = 1;
321 : call timab(607,1,tsec2)
322 : ABI_MALLOC(an_dev,(0:recpar%npt-1,0:nrec))
323 : ABI_MALLOC(bn2_dev,(0:recpar%npt-1,0:nrec))
324 : an_dev = zero
325 : bn2_dev = zero; bn2_dev(:,0) = one
326 :
327 : call cudarec( rset, exppot,an_dev,bn2_dev,&
328 : & beta,trotter,tolrec,dtset%recgratio,dtset%ngfft(:3),max_rec)
329 :
330 : max_rec = min(max_rec,nrec)
331 : a_wrk(0:max_rec,1:recpar%npt) = transpose(an_dev(0:,0:max_rec))
332 : b2_wrk(0:max_rec,1:recpar%npt) = transpose(bn2_dev(0:,0:max_rec))
333 :
334 : ABI_FREE(an_dev)
335 : ABI_FREE(bn2_dev)
336 : call timab(607,2,tsec2)
337 :
338 : ipointlocal = recpar%npt+1
339 :
340 : ! !DEBUG CUDA
341 : ! if(rset%debug)then
342 : ! if( rset%mpi%me==0)then
343 : ! do ipoint = 1,rset%par%npt,1
344 : ! kk=ipoint/(product(dtset%ngfft(:2)))
345 : ! jj=ipoint/dtset%ngfft(1)-kk*dtset%ngfft(2)
346 : ! ii=ipoint-jj*dtset%ngfft(1)-kk*dtset%ngfft(2)*dtset%ngfft(1)
347 : ! write(msg,'(a,4i8,2(a,a9,5f12.6))')&
348 : ! & 'pt',ipoint,ii,jj,kk,&
349 : ! & ch10,'an-gpu ',real(alocal(:4,ipoint)),&
350 : ! & ch10,'b2n-gpu ',real(b2local(:4,ipoint))
351 : ! call wrtout(std_out,msg,'COLL')
352 : ! end do
353 : ! endif
354 : ! endif
355 : ! !ENDDEBUG CUDA
356 : #endif
357 :
358 : else
359 38 : if (.not.(rset%tronc)) then
360 120 : graou1 : do kk = recpar%pt0%z,recpar%pt1%z,dtset%recgratio
361 1560 : do jj = 0,n2-1,dtset%recgratio
362 18778 : do ii = 0,n1-1,dtset%recgratio
363 17280 : ipoint = ii+(jj+kk*n2)*n1
364 : ! --Local position of atoms
365 17280 : if (ipoint<recpar%min_pt) cycle
366 : ! --Computation done by that proc
367 17280 : tim_fourdp=6
368 : call recursion(exppot,ii,jj,kk, &
369 : & a_wrk(:,ipointlocal), &
370 : & b2_wrk(:,ipointlocal), &
371 : & rho_wrk(ipointlocal),&
372 : & nrec, rset%efermi,tsmear,rtrotter,dim_trott, &
373 : & rset%ZT_p, &
374 : & tolrec,dtset%typat,rset%nl,&
375 : & rset%mpi,nfftrec,ngfftrec,rset%inf,&
376 17280 : & tim_fourdp,dtset%natom,projec,1)
377 17280 : ipointlocal = ipointlocal + 1
378 : ! write(std_out,*)'ipointlocal',ipoint,ipointlocal,ii,jj,kk
379 17280 : call prtwork(dtset%recgratio**3*ipointlocal*100*rset%mpi%nproc/nfftrec)
380 18689 : if(ipoint==recpar%max_pt) exit graou1
381 : end do
382 : end do
383 : end do graou1
384 : else !--We use a troncation
385 21 : ABI_MALLOC(exppotloc,(0:nfftrec-1))
386 21 : graou2 : do kk = recpar%pt0%z,recpar%pt1%z,dtset%recgratio
387 84 : do jj = 0,n2-1,dtset%recgratio
388 259 : do ii = 0,n1-1,dtset%recgratio
389 189 : ipoint = ii+(jj+kk*n2)*n1
390 189 : if (ipoint<recpar%min_pt) cycle
391 : ! computation done by that proc
392 12285 : exppotloc = zero
393 : ! --Traslation to move position on the ngfftrec grid center
394 756 : trasl = -(/ii,jj,kk/)+ngfftrec(:3)/2
395 189 : if(rset%nl%nlpsp) then
396 0 : do iatom=1,dtset%natom
397 : ! --local position of atoms
398 0 : gcart_loc(:,iatom) = rset%inf%gcart(:,iatom)+trasl
399 0 : gcart_loc(:,iatom) = modulo(gcart_loc(:,iatom),(/n1,n2,n3/))
400 : ! --Traslation of non-local projectors
401 0 : do ilmn = 1,rset%nl%lmnmax
402 0 : projec(:,:,:,ilmn,iatom) = reshape(rset%nl%projec(:,ilmn,dtset%typat(iatom)),shape=shape(projec(:,:,:,1,1)))
403 0 : do ii1=1,3
404 0 : projec(:,:,:,ilmn,iatom) = eoshift(projec(:,:,:,ilmn,iatom),shift=ngfftrec(ii1)/2-gcart_loc(ii1,iatom),dim=ii1)
405 : end do
406 : end do
407 : end do
408 : end if
409 :
410 : call reshape_pot(trasl,nfftf,nfftrec,&
411 : & dtset%ngfft(:3),ngfftrec(:3),&
412 189 : & exppot,exppotloc)
413 :
414 189 : tim_fourdp=6
415 : call recursion(exppotloc,ngfftrec(1)/2,ngfftrec(2)/2,ngfftrec(3)/2, &
416 : & a_wrk(:,ipointlocal), &
417 : & b2_wrk(:,ipointlocal), &
418 : & rho_wrk(ipointlocal),&
419 : & nrec, rset%efermi,tsmear,rtrotter,dim_trott, &
420 : & rset%ZT_p, &
421 : & tolrec,dtset%typat,rset%nl,&
422 : & rset%mpi,nfftrec,ngfftrec,rset%inf,&
423 189 : & tim_fourdp,dtset%natom,projec,1)
424 189 : ipointlocal = ipointlocal + 1
425 189 : call prtwork(factor*real(ipointlocal,dp))
426 245 : if(ipoint==recpar%max_pt) exit graou2
427 : end do
428 : end do
429 : end do graou2
430 7 : ABI_FREE(exppotloc)
431 : end if
432 :
433 : end if
434 38 : write(msg,'( a12,i12)')'ipointlocal',ipointlocal
435 38 : call wrtout(std_out,msg,'PERS')
436 38 : call timab(613+swt_tm,1,tsec2) !!--start time-counter: sync gpu-cpu
437 38 : call xmpi_barrier(rset%mpi%comm_bandfft)
438 38 : call timab(613+swt_tm,2,tsec2) !!--stop time-counter: sync gpu-cpu
439 :
440 : !!#############################################################
441 : !!--ASSIGNATION PARAMETERS TO MENAGE PARALLELISME OF TRANSGRID
442 38 : if((rset%load==1 .or. dtset%recgratio>1))then
443 : ! --Bufsize contains the values of the number of points calculated
444 : ! by any proc on the coarse grid; bufsize_f on the fine grid
445 7 : call timab(604,1,tsec2) !--start time-counter: transgrid
446 21 : ABI_MALLOC(bufsize,(0:rset%mpi%nproc-1))
447 14 : ABI_MALLOC(bufdispl,(0:rset%mpi%nproc-1))
448 14 : bufsize = 0;
449 7 : bufsize(rset%mpi%me) = rset%par%npt
450 7 : call xmpi_sum(bufsize,rset%mpi%comm_bandfft,ierr)
451 :
452 7 : bufdispl(0) = 0;
453 7 : if(rset%mpi%nproc>1) bufdispl(1:) = (/(sum(bufsize(:ii)),ii=0,rset%mpi%nproc-1)/)
454 14 : call timab(604,2,tsec2) !--stop time-counter: transgrid
455 : end if
456 : !!####################################################################
457 : !!--REDISTRIBUTION OF LOAD ON PROCS AFTER RECURSION IF CUDA IS USED
458 : #ifdef HAVE_GPU_CUDA
459 : if(rset%load==1)then
460 : call timab(604,1,tsec2) !--start time-counter: transgrid
461 : call xredistribute(rho_hyb,rset%GPU%par%vcount,rset%GPU%par%displs,&
462 : & rholocal,bufsize,bufdispl,rset%mpi%me,&
463 : & rset%mpi%nproc,&
464 : & rset%mpi%comm_bandfft,ierr)
465 :
466 :
467 : ABI_MALLOC(vcount_0,(0:rset%mpi%nproc-1))
468 : ABI_MALLOC(displs_0,(0:rset%mpi%nproc-1))
469 : ABI_MALLOC(vcount_1,(0:rset%mpi%nproc-1))
470 : ABI_MALLOC(displs_1,(0:rset%mpi%nproc-1))
471 :
472 : vcount_0 = 0
473 : vcount_0(rset%mpi%me) = rset%par%npt*(nrec+1)
474 : call xmpi_sum(vcount_0,rset%mpi%comm_bandfft,ierr)
475 : displs_0 = 0
476 : if(rset%mpi%nproc>1) displs_0(1:) = (/(sum(vcount_0(:ii)),ii=0,rset%mpi%nproc-1)/)
477 :
478 :
479 : vcount_1 = 0
480 : vcount_1(rset%mpi%me) = rset%GPU%par%npt*(nrec+1)
481 : call xmpi_sum(vcount_1,rset%mpi%comm_bandfft,ierr)
482 : displs_1 = 0
483 : if(rset%mpi%nproc>1) displs_1(1:) = (/(sum(vcount_1(:ii)),ii=0,rset%mpi%nproc-1)/)
484 :
485 :
486 : call xredistribute(a_hyb,vcount_1,displs_1,&
487 : & alocal,vcount_0,displs_0,&
488 : & rset%mpi%me,rset%mpi%nproc,&
489 : & rset%mpi%comm_bandfft,ierr)
490 :
491 : call xredistribute(b2_hyb,vcount_1,displs_1,&
492 : & b2local,vcount_0,displs_0,&
493 : & rset%mpi%me,rset%mpi%nproc,&
494 : & rset%mpi%comm_bandfft,ierr)
495 :
496 : nullify(rho_wrk,a_wrk,b2_wrk)
497 : ABI_FREE(rho_hyb)
498 : ABI_FREE(a_hyb)
499 : ABI_FREE(b2_hyb)
500 : ABI_FREE(vcount_0)
501 : ABI_FREE(displs_0)
502 : ABI_FREE(vcount_1)
503 : ABI_FREE(displs_1)
504 : call timab(604,2,tsec2) !--start time-counter: transgrid
505 : end if
506 : #endif
507 :
508 : !#############################################################
509 : !--TRANSGRID FOR THE DENSITY RHO AND THE COEFFICIENTS AN AND B2N
510 38 : if (dtset%recgratio>1) then
511 : ! --variables allocation and initialisation-------
512 7 : write (msg,'(a)')' - TRANSGRID USING -----'
513 7 : call wrtout(std_out,msg,'COLL')
514 7 : call timab(604,1,tsec2) !--start time-counter: transgrid
515 :
516 21 : ABI_MALLOC(rholocal_f,(rset%pawfgr%nfft))
517 21 : ABI_MALLOC(rhogf,(2,rset%pawfgr%nfft))
518 21 : ABI_MALLOC(rhogc,(2,rset%pawfgr%nfftc))
519 21 : ABI_MALLOC(rholoc_2,(1:rset%pawfgr%nfftc))
520 35 : ABI_MALLOC(ablocal_2,(1:rset%pawfgr%nfftc,0:nrec,2))
521 35 : ABI_MALLOC(ablocal_f,(1:rset%pawfgr%nfft,0:nrec,2))
522 35 : ABI_MALLOC(ablocal_1,(1:rset%par%npt,0:nrec,2))
523 :
524 7 : call rset%mpi%distribfft%free()
525 7 : call rset%mpi%distribfft%init('c',rset%mpi%nproc_fft,rset%pawfgr%ngfftc(2) ,rset%pawfgr%ngfftc(3))
526 7 : call rset%mpi%distribfft%init('f',rset%mpi%nproc_fft,rset%pawfgr%ngfft(2) ,rset%pawfgr%ngfft(3))
527 :
528 338933 : rholocal_f = zero; ablocal_f = zero; ablocal_2 = zero
529 :
530 2611 : ablocal_1(:,:,1) = transpose(alocal(:,1:rset%par%npt))
531 2611 : ablocal_1(:,:,2) = transpose(b2local(:,1:rset%par%npt))
532 :
533 7 : if(get_K_S_G==1 .and. dtset%recgratio>1 ) then
534 8 : ABI_MALLOC(aloc_copy,(0:nrec,1:rset%par%npt))
535 6 : ABI_MALLOC(b2loc_copy,(0:nrec,1:rset%par%npt))
536 436 : aloc_copy = alocal(:,1:rset%par%npt)
537 436 : b2loc_copy = b2local(:,1:rset%par%npt)
538 : end if
539 :
540 7 : if(rset%mpi%nproc ==1) then
541 : ! --SEQUENTIAL CASE--
542 203 : rholoc_2 = rholocal(1:rset%par%npt)
543 5236 : ablocal_2 = ablocal_1
544 : ! --Transigrid: coarse->fine
545 36862 : rhogf = zero; rhogc = zero
546 7 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,rset%pawfgr,rhogc,rhogf,rholoc_2,rholocal_f)
547 21 : do jj1 = 1,2
548 207 : do ipoint = 0,nrec
549 979476 : rhogf = zero; rhogc = zero
550 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,rset%pawfgr,rhogc,rhogf,&
551 200 : & ablocal_2(:,ipoint,jj1),ablocal_f(:,ipoint,jj1))
552 : end do
553 : end do
554 : ! --Assignation of the interpolated results on the fine grid--
555 12103 : rholocal = rholocal_f
556 172807 : alocal = transpose(ablocal_f(:,:,1))
557 172807 : b2local = transpose(ablocal_f(:,:,2))
558 :
559 : else
560 : ! --PARALLEL CASE--
561 0 : rholoc_2 = zero
562 : ! --Send on all procs rho,an,bn--
563 : call xmpi_allgatherv(rholocal(1:rset%par%npt),bufsize(rset%mpi%me),rholoc_2,&
564 0 : & bufsize,bufdispl,rset%mpi%comm_bandfft,ierr)
565 0 : do irec = 0,nrec
566 : call xmpi_allgatherv(ablocal_1(1:rset%par%npt,irec,1),bufsize(rset%mpi%me),&
567 0 : & ablocal_2(:,irec,1),bufsize,bufdispl,rset%mpi%comm_bandfft,ierr)
568 : call xmpi_allgatherv(ablocal_1(1:rset%par%npt,irec,2),bufsize(rset%mpi%me),&
569 0 : & ablocal_2(:,irec,2),bufsize,bufdispl,rset%mpi%comm_bandfft,ierr)
570 : end do
571 :
572 :
573 : ! --Transigrid: coarse->fine on differents procs (with respect
574 : ! the number of recursion)
575 0 : rhogf = zero; rhogc = zero
576 0 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,rset%pawfgr,rhogc,rhogf,rholoc_2,rholocal_f)
577 :
578 0 : do irec = 0,2*(nrec+1)-1
579 0 : ii1 = modulo(irec,rset%mpi%nproc)
580 0 : jj1 = 1+modulo(irec,2)
581 0 : kk1 = floor(irec/2.)
582 0 : if(maxval(abs(ablocal_2(:,kk1,jj1))) > tol10 .and. rset%mpi%me == ii1) then
583 0 : rhogf = zero; rhogc = zero
584 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,rset%pawfgr,rhogc,&
585 0 : & rhogf,ablocal_2(:,kk1,jj1),ablocal_f(:,kk1,jj1))
586 : end if
587 : end do
588 :
589 : ! --Recuperation of all interpolated results
590 : ! from procs to allprocs
591 0 : call xmpi_sum(ablocal_f,rset%mpi%comm_bandfft,ierr)
592 :
593 : ! --Assignation the interpolated results on the fine grid
594 : ! any procs to obtain the same point as in the standard recursion
595 0 : do ii1 = 0, rset%mpi%nproc-1
596 0 : jj1 = rset%par%displs(ii1)+1
597 0 : if(ii1 == rset%mpi%me) then
598 0 : alocal = transpose(ablocal_f(jj1:jj1+rset%par%vcount(ii1)-1,:,1))
599 0 : b2local = transpose(ablocal_f(jj1:jj1+rset%par%vcount(ii1)-1,:,2))
600 0 : rholocal = rholocal_f(jj1:jj1+rset%par%vcount(ii1)-1)
601 : end if
602 : end do
603 : end if
604 7 : ABI_FREE(ablocal_f)
605 7 : ABI_FREE(ablocal_2)
606 7 : ABI_FREE(ablocal_1)
607 7 : ABI_FREE(rhogf)
608 7 : ABI_FREE(rholocal_f)
609 7 : ABI_FREE(rhogc)
610 7 : ABI_FREE(rholoc_2)
611 :
612 7 : call timab(604,2,tsec2) !--stop time-counter: transgrid
613 : else
614 31 : write(msg,'(a)')' - TRANSGRID NOT USED --'
615 31 : call wrtout(std_out,msg,'COLL')
616 : end if
617 : !!--End transgrid
618 : !!###############################################################
619 : !###################################
620 : !--Fermi energy computation
621 : !--find the good mu by imposing the electrons number
622 : call fermisolverec(rset%efermi,rholocal,alocal,b2local,rset%debug,&
623 : & rset%min_nrec,tsmear,trotter,nelect,tol10,100, &
624 : & rset%par%ntranche,rset%mpi,inf_ucvol,& !.False. .and.&
625 76 : & (rset%tp==2 .or. rset%tp==3) .and. trotter>1)
626 :
627 : !#################################################################
628 : !######### ENTROPY AND GRAN POTENTIAL COMPUTATION ##################
629 38 : entropy = zero
630 38 : gran_pot = zero
631 38 : noentropie : if(get_K_S_G==1)then
632 11 : entropy1 = zero; entropy2 = zero ;entropy3 = zero; entropy4 = zero
633 11 : gran_pot1 = zero ; gran_pot2 = zero; gran_pot3 = zero; gran_pot4 = zero
634 :
635 : ! --Seek for the min of the path integral
636 11 : potmin = zero; nlpotmin = zero
637 17301 : if(dtset%rectesteg/=1) potmin = minval(vtrial(:,1))
638 11 : if(rset%nl%nlpsp) nlpotmin = minval(rset%nl%eival(:,:,:))
639 11 : xmax = exp(-ratio2*(potmin+nlpotmin-rset%efermi))
640 :
641 11 : dim_entro = 0; if(rset%debug) dim_entro = 4;
642 :
643 11 : if(dtset%recgratio>1) then
644 : ! --Recgratio>1
645 6 : ABI_MALLOC(rhogf,(2,rset%pawfgr%nfft))
646 6 : ABI_MALLOC(rhogc,(2,rset%pawfgr%nfftc))
647 6 : ABI_MALLOC(entropy_v_f,(rset%pawfgr%nfft,0:4))
648 6 : ABI_MALLOC(entropy_v_c,(rset%pawfgr%nfftc,0:4))
649 6 : ABI_MALLOC(entropy_v_2,(1:rset%par%npt,0:4))
650 4 : ABI_MALLOC(gran_pot_v_f,(rset%pawfgr%nfft,0:4))
651 4 : ABI_MALLOC(gran_pot_v_c,(rset%pawfgr%nfftc,0:4))
652 4 : ABI_MALLOC(gran_pot_v_2,(1:rset%par%npt,0:4))
653 :
654 17852 : entropy_v_c = zero; entropy_v_f = zero; entropy_v_2 = zero
655 17852 : gran_pot_v_c = zero; gran_pot_v_f = zero; gran_pot_v_2 = zero
656 :
657 :
658 56 : do ipoint = 1,rset%par%npt
659 : call entropyrec(exp(rset%efermi*ratio2)*aloc_copy(:,ipoint), &
660 : & exp(rset%efermi*ratio1)*b2loc_copy(:,ipoint), &
661 : & nrec,trotter,entropy_v_2(ipoint,0),two,&
662 : & rset%debug,n_pt_integ_entropy,perc_vmin*xmax,&
663 : & entropy_v_2(ipoint,1),&
664 : & entropy_v_2(ipoint,2),&
665 : & entropy_v_2(ipoint,3),&
666 810 : & entropy_v_2(ipoint,4))
667 :
668 : call gran_potrec(exp(rset%efermi*ratio2)*aloc_copy(:,ipoint), &
669 : & exp(rset%efermi*ratio1)*b2loc_copy(:,ipoint), &
670 : & nrec,trotter,gran_pot_v_2(ipoint,0),two,&
671 : & rset%debug,n_pt_integ_entropy,perc_vmin*xmax,&
672 : & gran_pot_v_2(ipoint,1),&
673 : & gran_pot_v_2(ipoint,2),&
674 : & gran_pot_v_2(ipoint,3),&
675 812 : & gran_pot_v_2(ipoint,4))
676 : end do
677 2 : ABI_FREE(aloc_copy)
678 2 : ABI_FREE(b2loc_copy)
679 :
680 2 : call timab(613+swt_tm,1,tsec2) !!--start time-counter: sync gpu-cpu
681 2 : call xmpi_barrier(rset%mpi%comm_bandfft)
682 2 : call timab(613+swt_tm,2,tsec2) !!--stop time-counter: sync gpu-cpu
683 :
684 2 : call timab(604,1,tsec2) !--start time-counter: transgrid
685 2 : if(rset%mpi%nproc==1) then
686 284 : entropy_v_c = entropy_v_2
687 284 : gran_pot_v_c = gran_pot_v_2
688 : end if
689 4 : do ii1 = 0,dim_entro
690 : call xmpi_allgatherv(entropy_v_2(:,ii1),bufsize(rset%mpi%me),&
691 : & entropy_v_c(:,ii1),bufsize,bufdispl,&
692 2 : & rset%mpi%comm_bandfft,ierr)
693 : call xmpi_allgatherv(gran_pot_v_2(:,ii1),bufsize(rset%mpi%me),&
694 : & gran_pot_v_c(:,ii1),bufsize,bufdispl,&
695 2 : & rset%mpi%comm_bandfft,ierr)
696 :
697 58 : if(maxval(abs(entropy_v_c(:,ii1))) > tol10) then
698 10532 : rhogf = zero; rhogc = zero
699 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,&
700 2 : & rset%pawfgr,rhogc,rhogf,entropy_v_c(:,ii1),entropy_v_f(:,ii1))
701 : end if
702 60 : if(maxval(abs(gran_pot_v_c(:,ii1))) >tol10) then
703 10532 : rhogf = zero; rhogc = zero
704 : call transgrid(1,rset%mpi,dtset%nspden,1,0,0,1,&
705 2 : & rset%pawfgr,rhogc,rhogf,gran_pot_v_c(:,ii1),gran_pot_v_f(:,ii1))
706 : end if
707 : end do
708 2 : call timab(604,2,tsec2) !--stop time-counter: transgrid
709 :
710 3458 : entropy = sum(entropy_v_f(:,0))
711 3458 : gran_pot = sum(gran_pot_v_f(:,0))
712 :
713 2 : if(rset%debug)then
714 0 : entropy1 = sum(entropy_v_f(:,1))
715 0 : entropy2 = sum(entropy_v_f(:,2))
716 0 : entropy3 = sum(entropy_v_f(:,3))
717 0 : entropy4 = sum(entropy_v_f(:,4))
718 0 : gran_pot1 = sum(gran_pot_v_f(:,1))
719 0 : gran_pot2 = sum(gran_pot_v_f(:,2))
720 0 : gran_pot3 = sum(gran_pot_v_f(:,3))
721 0 : gran_pot4 = sum(gran_pot_v_f(:,4))
722 : end if
723 :
724 2 : ABI_FREE(entropy_v_f)
725 2 : ABI_FREE(entropy_v_c)
726 2 : ABI_FREE(entropy_v_2)
727 2 : ABI_FREE(rhogf)
728 2 : ABI_FREE(rhogc)
729 2 : ABI_FREE(gran_pot_v_f)
730 2 : ABI_FREE(gran_pot_v_c)
731 2 : ABI_FREE(gran_pot_v_2)
732 :
733 :
734 : else
735 : ! --Recgratio=1
736 5193 : do ipoint = 1,rset%par%ntranche
737 : call entropyrec(exp(rset%efermi*ratio2)*alocal(:,ipoint), &
738 : & exp(rset%efermi*ratio1)*b2local(:,ipoint), &
739 : & nrec,trotter,entropylocal,two,&
740 : & rset%debug,n_pt_integ_entropy,perc_vmin*xmax,&
741 : & entropylocal1,entropylocal2,&
742 77760 : & entropylocal3,entropylocal4)
743 : call gran_potrec(exp(rset%efermi*ratio2)*alocal(:,ipoint), &
744 : & exp(rset%efermi*ratio1)*b2local(:,ipoint), &
745 : & nrec,trotter,gran_pot_local,two,& !/ucvol,&
746 : & rset%debug,n_pt_integ_entropy,perc_vmin*xmax,&
747 : & gran_pot_local1,gran_pot_local2,&
748 77760 : & gran_pot_local3,gran_pot_local4)
749 :
750 5184 : entropy = entropy + entropylocal
751 5184 : gran_pot = gran_pot + gran_pot_local
752 5193 : if(rset%debug)then
753 0 : entropy1 = entropy1 + entropylocal1
754 0 : entropy2 = entropy2 + entropylocal2
755 0 : entropy3 = entropy3 + entropylocal3
756 0 : entropy4 = entropy4 + entropylocal4
757 0 : gran_pot1 = gran_pot1 + gran_pot_local1
758 0 : gran_pot2 = gran_pot2 + gran_pot_local2
759 0 : gran_pot3 = gran_pot3 + gran_pot_local3
760 0 : gran_pot4 = gran_pot4 + gran_pot_local4
761 : end if
762 :
763 : end do
764 :
765 9 : call xmpi_sum(entropy,rset%mpi%comm_bandfft ,ierr)
766 9 : call xmpi_sum(gran_pot,rset%mpi%comm_bandfft ,ierr)
767 9 : if(rset%debug)then
768 0 : call xmpi_sum(entropy1,rset%mpi%comm_bandfft ,ierr)
769 0 : call xmpi_sum(entropy2,rset%mpi%comm_bandfft ,ierr)
770 0 : call xmpi_sum(entropy3,rset%mpi%comm_bandfft ,ierr)
771 0 : call xmpi_sum(entropy4,rset%mpi%comm_bandfft ,ierr)
772 0 : call xmpi_sum(gran_pot1,rset%mpi%comm_bandfft ,ierr)
773 0 : call xmpi_sum(gran_pot2,rset%mpi%comm_bandfft ,ierr)
774 0 : call xmpi_sum(gran_pot3,rset%mpi%comm_bandfft ,ierr)
775 0 : call xmpi_sum(gran_pot4,rset%mpi%comm_bandfft ,ierr)
776 : end if
777 : end if
778 :
779 11 : if(rset%debug)then
780 : write(msg,'(2(2a,4(2a,es11.4,a)))')&
781 0 : & ' --------------------------' ,ch10, &
782 0 : & ' entropy, horiz path=',' ',entropy1,ch10, &
783 0 : & ' entropy, xmax path=',' ',entropy2,ch10, &
784 0 : & ' entropy, xmin path=',' ',entropy3,ch10, &
785 0 : & ' entropy, zero path=',' ',entropy4,ch10, &
786 0 : & ' --------------------------' ,ch10, &
787 0 : & ' -omega/T, horiz path=',' ',gran_pot1,ch10, &
788 0 : & ' -omega/T, xmax path=',' ',gran_pot2,ch10, &
789 0 : & ' -omega/T, xmin path=',' ',gran_pot3,ch10, &
790 0 : & ' -omega/T, zero path=',' ',gran_pot4,ch10
791 0 : call wrtout(std_out,msg,'COLL')
792 : end if
793 :
794 11 : e_eigenvalues=tsmear*(entropy-gran_pot) + rset%efermi*nelect
795 : ! --In reality gran_pot is not the gran potential but the
796 : ! potential omega=-PV (Landau-potential or grand-potential)
797 : ! divided by -T so the internal energy
798 : ! U:=e_eigenvalues= TS+omega+muN = ST-T*sum(ln(1-n))+muN =
799 : ! T(S-gran_pot)+muN
800 :
801 :
802 11 : if(rset%nl%nlpsp) then
803 : call nlenergyrec(rset,enlx,exppot,dtset%ngfft,dtset%natom,&
804 0 : & dtset%typat,tsmear,trotter,tolrec)
805 : end if
806 : end if noentropie
807 : !##### END ENTROPY AND GRAN POTENTIAL COMPUTATION ##################
808 : !#################################################################
809 :
810 : !if(associated(projec))
811 38 : if(rset%nl%nlpsp) then
812 0 : ABI_FREE(projec)
813 : end if
814 38 : if(associated(gcart_loc)) then
815 35 : ABI_FREE(gcart_loc)
816 : end if
817 38 : if((dtset%recgratio/=1 .or. rset%load==1)) then
818 7 : ABI_FREE(bufdispl)
819 7 : ABI_FREE(bufsize)
820 : end if
821 : !------------------------------------------------------------------
822 : !--Check if the convergence is reached for rho
823 29452 : drho = maxval(abs(rhor(min_pt:max_pt,1)-rholocal(:)))
824 : drhomax = drho
825 38 : call xmpi_max(drho,drhomax,rset%mpi%comm_bandfft,ierr)
826 :
827 : !write(std_out,*)'drhomax,toldrho',drhomax,toldrho
828 38 : if(drhomax<toldrho)then
829 22 : rset%quitrec = rset%quitrec+1
830 : else
831 16 : rset%quitrec = 0
832 : end if
833 :
834 : !-------------------------------------------------------------------
835 : !--Density on all procs
836 29414 : rhor(min_pt:max_pt,1) = rholocal(:)
837 38 : if(rset%mpi%nproc /= 1)then
838 : call xmpi_allgatherv(rholocal,rset%par%ntranche,rhor(:,1),&
839 : & rset%par%vcount,rset%par%displs,&
840 28 : & rset%mpi%comm_band,ierr)
841 : end if
842 :
843 : !--------------------------------------------------------------------
844 : !--2nd EKIN CALCULATION: this method is used
845 38 : noekin2 : if(get_K_S_G==1)then
846 8651 : intrhov = (inf_ucvol)*sum(rholocal*vtrial(min_pt:max_pt,1))
847 11 : call xmpi_sum(intrhov,rset%mpi%comm_bandfft ,ierr)
848 :
849 11 : ek = e_eigenvalues-intrhov-enlx
850 :
851 :
852 11 : if(rset%debug) then
853 0 : write (msg,'(2a,3f15.10,2a,3f15.10,2a,f15.10,a)') ch10,&
854 0 : & ' ek,int(rho*V),ek+int(rho*V) ', ek, intrhov, ek+ intrhov,ch10, &
855 0 : & ' kT*S, kT*sum(ln(...)), diff ', tsmear*entropy, tsmear*gran_pot, tsmear*(entropy-gran_pot),ch10, &
856 0 : & ' kT(S-sum(ln(...)))+mu*nelect', tsmear*(entropy-gran_pot)+rset%efermi*nelect,ch10
857 0 : call wrtout(std_out,msg,'COLL')
858 : end if
859 :
860 :
861 : end if noekin2
862 : !--------------------------------------------------------------------
863 38 : fermie = rset%efermi
864 :
865 : !--------------------------------------------------------
866 : !!--At the first step to find the max number of recursion
867 : !! needed to convergence, then redefine nrec.
868 38 : if(initialized==0 .and. dtset%ntime>0) then
869 6 : call Calcnrec(rset,b2local)
870 : end if
871 :
872 : !--------------------------------------------------------
873 38 : call rset%mpi%distribfft%free()
874 38 : call rset%mpi%distribfft%init('c',rset%mpi%nproc_fft,rset%ngfftrec(2),rset%ngfftrec(3))
875 38 : call rset%mpi%distribfft%init('f',rset%mpi%nproc_fft,dtset%ngfft(2),dtset%ngfft(3))
876 :
877 : !--Printing results
878 : write(msg,'(3a,f15.10)')&
879 38 : & ' -- Results: --------------------------------------',ch10,&
880 76 : & ' mu =',rset%efermi
881 38 : call wrtout(std_out,msg,'COLL')
882 38 : if(get_K_S_G==1)then
883 : write(msg,'(a,f15.10,6(2a,f15.10))')&
884 11 : & ' potmin =',potmin,ch10,&
885 11 : & ' <V_eff> =',intrhov,ch10,&
886 11 : & ' entropy =',entropy,ch10,&
887 11 : & ' -omega/T =',gran_pot,ch10,&
888 11 : & ' eigenvalues =',e_eigenvalues,ch10,&
889 11 : & ' kinetic =',ek,ch10,&
890 22 : & ' non-loc ene =',enlx
891 11 : call wrtout(std_out,msg,'COLL')
892 : end if
893 1938 : write(msg,'(a,50a)')' ',('-',ii=1,50)
894 38 : call wrtout(std_out,msg,'COLL')
895 : !write(std_out,*)'is the pressure ',gran_pot*tsmear/(rset%inf%ucvol*real(nfftrec,dp))
896 :
897 : !--Structured debugging : if rset%debug=T, stop here.
898 : if(.false.)then !(rset%debug)
899 : call wrtout(std_out,' rhor ','PERS')
900 : write(std_out,*)rhor(:,1)
901 : call wrtout(std_out,' ','COLL')
902 : write(msg,'(a,2d10.3)')' temps recursion ',tsec
903 : call wrtout(std_out,msg,'COLL')
904 : write(msg,'(a,l1,a)') ' vtorhorec : rset%debug=-',rset%debug,', debugging mode => stop '
905 : ABI_ERROR(msg)
906 : end if
907 :
908 38 : call timab(600,2,tsec2)
909 38 : call timab(21,2,tsec)
910 :
911 : call symrhg(1,gprimd,irrzon,rset%mpi,nfftf,&
912 : & dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),dtset%ngfft,dtset%nspden,&
913 38 : & dtset%nsppol,dtset%nsym,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
914 :
915 76 : end subroutine vtorhorec
916 : !!***
917 :
918 : !!****f* ABINIT/entropyrec
919 : !! NAME
920 : !! entropyrec
921 : !!
922 : !! FUNCTION
923 : !! This routine computes the local part of the entropy at a point using a path integral,
924 : !! in the recursion method.
925 : !!
926 : !! an, bn2 : coefficient given by the recursion.
927 : !! nrec=order of recursion
928 : !! trotter=trotter parameter
929 : !! multce=a multiplicator for computing entropy ; 2 for non-spin-polarized system
930 : !! debug_rec=debug variable
931 : !! n_pt_integ=number of points of integration for the path integral
932 : !! xmax =max point of integration on the real axis
933 :
934 : !! OUTPUT
935 : !! ent_out=entropy at the point
936 : !! ent_out1,ent_out2,ent_out3,ent_out4=debug entropy at the point
937 : !!
938 : !! NOTES
939 : !! at this time :
940 : !! - multce should be not used
941 : !! - the routine should be integraly rewrited and use the routine recursion.
942 : !! - only modified for p /= 0
943 : !!
944 : !! SOURCE
945 :
946 5238 : subroutine entropyrec(an,bn2,nrec,trotter,ent_out,multce,debug_rec, &
947 : & n_pt_integ,xmax,&
948 : & ent_out1,ent_out2,ent_out3,ent_out4)
949 :
950 : !Arguments -------------------------------
951 : !scalars
952 : integer,intent(in) :: n_pt_integ,nrec,trotter
953 : logical,intent(in) :: debug_rec
954 : real(dp), intent(in) :: multce,xmax
955 : real(dp),intent(out) :: ent_out,ent_out1,ent_out2,ent_out3,ent_out4
956 : !arrays
957 : real(dp),intent(in) :: an(0:nrec),bn2(0:nrec)
958 :
959 : !Local variables-------------------------------
960 : !scalars
961 : integer, parameter :: level = 7
962 : integer, save :: first_en = 1
963 : integer :: ii,kk,n_pt_integ_path2,n_pt_integ_path3
964 : real(dp) :: arg,epsilo,step,twotrotter,xmin,dr_step
965 : complex(dp) :: D,Dnew,Dold,N,Nnew,Nold,dz_path,ent_acc,ent_acc1,ent_acc2
966 : complex(dp) :: ent_acc3,ent_acc4
967 : complex(dp) :: funczero,z_path,zj
968 : complex(dp) ::delta_calc
969 : character(len=500) :: msg
970 : !arrays
971 : real(dp) :: tsec(2)
972 : real(dp) :: iif,factor
973 : ! *************************************************************************
974 :
975 5238 : call timab(610,1,tsec)
976 :
977 : !structured debugging if debug_rec=T : print detailled result the first time we enter entropyrec
978 :
979 5238 : if(debug_rec .and. first_en==1)then
980 0 : write(msg,'(a)')' '
981 0 : call wrtout(std_out,msg,'PERS')
982 0 : write(msg,'(a)')' entropyrec : enter '
983 0 : call wrtout(std_out,msg,'PERS')
984 0 : write(msg,'(a,i6)')'n_pt_integ ' , n_pt_integ
985 0 : call wrtout(std_out,msg,'COLL')
986 : end if
987 :
988 5238 : ent_out = zero
989 5238 : ent_out1 = zero
990 5238 : ent_out2 = zero
991 5238 : ent_out3 = zero
992 5238 : ent_out4 = zero
993 5238 : ent_acc = czero
994 5238 : ent_acc1 = czero
995 5238 : ent_acc2 = czero
996 5238 : ent_acc3 = czero
997 5238 : ent_acc4 = czero
998 :
999 : !path parameters
1000 5238 : twotrotter = max(two*real(trotter,dp),one)
1001 5238 : if(trotter==0)then
1002 : factor = tol5
1003 : arg =pi*three_quarters
1004 : zj = cmplx(-one,one-sin(arg),dp)
1005 : else
1006 5238 : factor = xmax/ten
1007 5238 : arg = pi/twotrotter
1008 5238 : zj = cmplx( cos(arg) , sin(arg),dp )
1009 : end if
1010 :
1011 5238 : epsilo = factor*sin( arg )
1012 5238 : xmin = factor*cos( arg )
1013 5238 : step = (xmax-xmin)/real(n_pt_integ,dp)
1014 :
1015 : !####################################################################
1016 : ![xmax + i*epsilo,xmin + i*epsilo]
1017 5238 : dr_step = one/real(n_pt_integ,dp)
1018 534276 : path1: do ii = 0,n_pt_integ
1019 529038 : z_path = cmplx(xmin+real(ii,dp)*(xmax-xmin)*dr_step,epsilo,dp)
1020 529038 : dz_path = -cmplx((xmax-xmin)*dr_step,zero,dp)
1021 :
1022 529038 : Nold = czero
1023 529038 : Dold = cone
1024 529038 : N = cone
1025 529038 : D = z_path - cmplx(an(0),zero,dp)
1026 :
1027 2116152 : do kk=1,nrec
1028 2116152 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1029 2116152 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1030 :
1031 2116152 : Nold = N
1032 2116152 : Dold = D
1033 2116152 : N = Nnew
1034 2116152 : D = Dnew
1035 :
1036 2645190 : if(kk/=nrec)then
1037 2116152 : if((bn2(kk+1)<tol14))exit
1038 : end if
1039 : end do
1040 :
1041 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1042 529038 : delta_calc = func1_rec(z_path**twotrotter)*(N/D)*dz_path
1043 534276 : if(ii==0.or.ii==n_pt_integ)then
1044 10476 : ent_acc = ent_acc + half*delta_calc
1045 10476 : ent_acc1 = ent_acc1 + half*delta_calc
1046 : else
1047 518562 : ent_acc = ent_acc + delta_calc
1048 518562 : ent_acc1 = ent_acc1 + delta_calc
1049 : end if
1050 : end do path1
1051 :
1052 :
1053 : !####################################################################
1054 : ![1/2zj,0]
1055 5238 : if(epsilo/step>100.d0)then
1056 0 : n_pt_integ_path2 = int((factor*abs(zj))/step)+1
1057 : else
1058 5238 : n_pt_integ_path2 = 100
1059 : end if
1060 :
1061 5238 : if(trotter/=0)then
1062 5238 : n_pt_integ_path3 = 0
1063 5238 : dr_step = one/real(n_pt_integ_path2,dp)
1064 5238 : dz_path = -cmplx(xmin,epsilo,dp)*dr_step
1065 534276 : path5: do ii = 0,n_pt_integ_path2
1066 529038 : z_path = cmplx(real(ii,dp)*xmin,real(ii,dp)*epsilo,dp)*dr_step
1067 534276 : if(abs(z_path)>tol14)then
1068 523800 : Nold = czero
1069 523800 : Dold = cone
1070 523800 : N = cone
1071 523800 : D = z_path - cmplx(an(0),zero,dp)
1072 2095200 : do kk=1,nrec
1073 2095200 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1074 2095200 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1075 2095200 : Nold = N
1076 2095200 : Dold = D
1077 2095200 : N = Nnew
1078 2095200 : D = Dnew
1079 2619000 : if(kk/=nrec)then
1080 2095200 : if((bn2(kk+1)<tol14))exit
1081 : end if
1082 : end do
1083 :
1084 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1085 523800 : if(abs(z_path)**twotrotter>tiny(one)) then
1086 523800 : funczero = func1_rec(z_path**twotrotter)
1087 : else
1088 : funczero = czero
1089 : end if
1090 523800 : delta_calc = funczero*N/D*dz_path
1091 523800 : if(ii==0.or.ii==n_pt_integ_path2)then
1092 5238 : ent_acc = ent_acc + half*delta_calc
1093 5238 : if(debug_rec) ent_acc3 = ent_acc3 + half*delta_calc
1094 : else
1095 518562 : ent_acc = ent_acc + funczero*delta_calc
1096 518562 : if(debug_rec) ent_acc3 = ent_acc3 + funczero*delta_calc
1097 : end if
1098 : end if
1099 : end do path5
1100 :
1101 : else ! trotter==0
1102 :
1103 0 : n_pt_integ_path3 = max(100,int((epsilo*half*pi)/real(step,dp))+1)
1104 0 : dr_step = one/real(n_pt_integ_path3,dp)
1105 0 : path6: do ii = 0,n_pt_integ_path3
1106 0 : iif=half*pi*real(ii,dp)*dr_step
1107 0 : z_path = epsilo*cmplx(-cos(iif),1-sin(iif),dp)
1108 0 : dz_path = epsilo*cmplx(sin(iif),-cos(iif),dp)*half*pi*dr_step
1109 0 : if(abs(z_path)**twotrotter>tol14)then
1110 0 : Nold = czero
1111 0 : Dold = cone
1112 0 : N = cone
1113 0 : D = z_path - cmplx(an(0),zero,dp)
1114 0 : do kk=1,nrec
1115 0 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1116 0 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1117 0 : Nold = N
1118 0 : Dold = D
1119 0 : N = Nnew
1120 0 : D = Dnew
1121 0 : if(kk/=nrec .and. bn2(kk+1)<tol14) exit !-EXIT
1122 : end do
1123 :
1124 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1125 0 : delta_calc = func1_rec(z_path**twotrotter) * N/D * dz_path
1126 0 : if(ii==0.or.ii==n_pt_integ_path3)then
1127 0 : ent_acc = ent_acc + half*delta_calc
1128 0 : if(debug_rec) ent_acc3 = ent_acc3 + half*delta_calc
1129 : else
1130 0 : ent_acc = ent_acc + delta_calc !<r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1131 0 : if(debug_rec) ent_acc3 = ent_acc3 + delta_calc !<r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1132 : end if
1133 : end if
1134 : end do path6
1135 :
1136 : end if
1137 :
1138 5238 : if(first_en==1 .and. debug_rec) then
1139 : write(msg,'(a,i5,2a,i5,2a,i5,2a,es11.4,2a,es11.4,2a,es11.4)')&
1140 0 : & 'n_pt_path =',n_pt_integ,ch10,&
1141 0 : & 'n_pt_path2 =',n_pt_integ_path2,ch10,&
1142 0 : & 'n_pt_path3 =',n_pt_integ_path3,ch10,&
1143 0 : & 'xmin =',xmin,ch10,&
1144 0 : & 'xmax =',xmax,ch10,&
1145 0 : & 'epsilon =',epsilo
1146 0 : call wrtout(std_out,msg,'COLL')
1147 0 : first_en = 0
1148 : end if
1149 :
1150 : !####################################################################
1151 : ![xmax,xmax+i*epsilo]
1152 5238 : dr_step = one/real(n_pt_integ_path2,dp)
1153 5238 : dz_path = cmplx(zero,epsilo*dr_step,dp)
1154 534276 : path4: do ii = 0,n_pt_integ_path2
1155 529038 : z_path = cmplx(xmax,real(ii,dp)*epsilo*dr_step,dp)
1156 :
1157 529038 : Nold = czero
1158 529038 : Dold = cone
1159 529038 : N = cone
1160 529038 : D = z_path - cmplx(an(0),zero,dp)
1161 :
1162 2116152 : do kk=1,nrec
1163 2116152 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1164 2116152 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1165 :
1166 2116152 : Nold = N
1167 2116152 : Dold = D
1168 2116152 : N = Nnew
1169 2116152 : D = Dnew
1170 :
1171 2645190 : if(kk/=nrec)then
1172 2116152 : if((bn2(kk+1)<tol14))exit
1173 : end if
1174 : end do
1175 :
1176 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1177 529038 : delta_calc = func1_rec(z_path**twotrotter)*N/D*dz_path
1178 534276 : if(ii==0.or.ii==n_pt_integ_path2)then
1179 :
1180 10476 : ent_acc = ent_acc + half*delta_calc
1181 10476 : if(debug_rec) ent_acc2 = ent_acc2 + half*delta_calc
1182 : else
1183 518562 : ent_acc = ent_acc + delta_calc
1184 518562 : if(debug_rec) ent_acc2 = ent_acc2 + delta_calc
1185 : end if
1186 : end do path4
1187 :
1188 :
1189 5238 : ent_out = multce*real(ent_acc*cmplx(zero,-piinv,dp),dp)
1190 5238 : if(debug_rec) then
1191 0 : ent_out1 = multce*real(ent_acc1*cmplx(zero,-piinv,dp),dp)
1192 0 : ent_out2 = multce*real(ent_acc2*cmplx(zero,-piinv,dp),dp)
1193 0 : ent_out3 = multce*real(ent_acc3*cmplx(zero,-piinv,dp),dp)
1194 0 : ent_out4 = multce*real(ent_acc4*cmplx(zero,-piinv,dp),dp)
1195 : end if
1196 :
1197 5238 : call timab(610,2,tsec)
1198 :
1199 : contains
1200 :
1201 : !function to integrate over the path
1202 : !func1_rec(z_path,twotrotter) = ( z_path**twotrotter/(1+z_path**twotrotter)*log(1+1/z_path**twotrotter)+& !- f*ln(f)
1203 : !&1/(1+z_path**twotrotter)*log(1+z_path**twotrotter)) !- (1-f)*ln(1-f)
1204 :
1205 : !func1_rec(z_path_pow) = z_path_pow/(cone+z_path_pow)*log(cone+cone/z_path_pow)+& !- f*ln(f)
1206 : !&cone/(cone+z_path_pow)*log(cone+z_path_pow) !- (1-f)*ln(1-f)
1207 :
1208 : !other expression of func for a path like ro(t)*exp(2*i*pi/(2*p)*(j+1/2))
1209 :
1210 1581876 : function func1_rec(z)
1211 :
1212 : complex(dp) :: func1_rec
1213 : complex(dp),intent(in) :: z
1214 :
1215 1581876 : func1_rec = z/(cone+z)*log(cone+cone/z)+ cone/(cone+z)*log(cone+z)
1216 :
1217 1581876 : end function func1_rec
1218 :
1219 : end subroutine entropyrec
1220 : !!***
1221 :
1222 : !!****f* ABINIT/fermisolverec
1223 : !! NAME
1224 : !! fermisolverec
1225 : !!
1226 : !! FUNCTION
1227 : !! This routine computes the fermi energy in order to have a given number of
1228 : !! valence electrons in the recursion method, using a Ridder s Method
1229 : !!
1230 : !! INPUTS
1231 : !! debug_rec=debugging variable
1232 : !! nb_rec=order of recursion
1233 : !! nb_point=number of discretization point in one dimension (=n1=n2=n3)
1234 : !! temperature=temperature (Hartree)
1235 : !! trotter=trotter parameter
1236 : !! nelect=number of valence electrons (dtset%nelect)
1237 : !! acc=accuracy for the fermi energy
1238 : !! max_it=maximum number of iteration for the Ridder's Method
1239 : !! long_tranche=number of point computed by thi proc
1240 : !! mpi_enreg=information about MPI parallelization
1241 : !! inf_ucvol=infinitesimal unit cell volume
1242 : !! gputopo=true if topology gpu-cpu= 2 or 3
1243 : !!
1244 : !! OUTPUT
1245 : !!
1246 : !! SIDE EFFECTS
1247 : !! fermie=fermi energy
1248 : !! rho=density, recomputed for the new fermi energy
1249 : !! a, b2 : coefficient given by recursion recomputed for the new fermi energy
1250 : !!
1251 : !! NOTES
1252 : !! at this time :
1253 : !!
1254 : !! SOURCE
1255 :
1256 38 : subroutine fermisolverec(fermie,rho,a,b2,debug_rec,nb_rec, &
1257 : & temperature,trotter,nelect, &
1258 : & acc, max_it, &
1259 : & long_tranche,mpi_enreg,&
1260 : & inf_ucvol,gputopo)
1261 :
1262 : !Arguments -------------------------------
1263 : !scalars
1264 : integer,intent(in) :: long_tranche,max_it,nb_rec,trotter
1265 : logical,intent(in) :: debug_rec,gputopo
1266 : real(dp),intent(in) :: acc,inf_ucvol,nelect,temperature
1267 : real(dp), intent(inout) :: fermie
1268 : type(MPI_type),intent(in) :: mpi_enreg
1269 : !arrays
1270 : real(dp), intent(inout) :: a(0:nb_rec,long_tranche), b2(0:nb_rec,long_tranche)
1271 : real(dp), intent(inout) :: rho(long_tranche)
1272 :
1273 : !Local variables-------------------------------
1274 : !scalars
1275 : integer :: ierr,ii,ipointlocal,nn,dim_trott
1276 : real(dp) :: beta,fermieh,fermiel,fermiem,fermienew,nelecth,nelectl,nelectm
1277 : real(dp) :: nelectnew,res_nelecth,res_nelectl,res_nelectm,res_nelectnew
1278 : real(dp) :: rtrotter,ss,fermitol
1279 : character(len=500) :: msg
1280 : !arrays
1281 : real(dp) :: tsec(2)
1282 76 : real(dp) :: rhotry(long_tranche)
1283 : !no_abirules
1284 : #ifdef HAVE_GPU_CUDA
1285 : integer :: swt_tm,npitch
1286 : real(cudap) :: rhocu(long_tranche)
1287 : real(dp) :: tsec2(2)
1288 : #endif
1289 :
1290 : ! *************************************************************************
1291 :
1292 : #ifdef HAVE_GPU_CUDA
1293 : swt_tm = 0
1294 : #endif
1295 :
1296 38 : call timab(609,1,tsec)
1297 :
1298 38 : beta = one/temperature
1299 38 : rtrotter = max(half,real(trotter,dp))
1300 38 : dim_trott = max(0,2*trotter-1)
1301 :
1302 38 : write(msg,'(a)')' -- fermisolverec ---------------------------------'
1303 38 : call wrtout(std_out,msg,'COLL')
1304 38 : if(debug_rec) then
1305 0 : write (msg,'(a,d10.3)')' nelect= ',nelect
1306 0 : call wrtout(std_out,msg,'COLL')
1307 : end if
1308 : !initialisation of fermiel
1309 38 : fermiel = fermie
1310 38 : call timab(609,2,tsec)
1311 :
1312 : !initialisation fermitol
1313 38 : fermitol = acc
1314 : #ifdef HAVE_GPU_CUDA_SP
1315 : if(gputopo) fermitol = 1.d-3
1316 : #endif
1317 :
1318 38 : if(gputopo) then
1319 : #ifdef HAVE_GPU_CUDA
1320 : swt_tm = 1
1321 : ! allocate array an and bn2 on gpu for computation of trotter formula
1322 : call alloc_dens_cuda(long_tranche,nb_rec,dim_trott,npitch,&
1323 : & real(a,cudap),real(b2,cudap))
1324 :
1325 : call timab(617,1,tsec)
1326 : call density_cuda(npitch,long_tranche,nb_rec,dim_trott,&
1327 : & real(fermiel,cudap),real(temperature,cudap),&
1328 : & real(rtrotter,cudap),real(inf_ucvol,cudap),&
1329 : & real(tol14,cudap),&
1330 : & rhocu)
1331 : rhotry = real(rhocu,dp)
1332 : call timab(617,2,tsec)
1333 : #endif
1334 : else
1335 29414 : do ipointlocal = 1,long_tranche
1336 : call density_rec(a(:,ipointlocal),&
1337 : & b2(:,ipointlocal),&
1338 : & rhotry(ipointlocal),&
1339 : & nb_rec,fermiel,temperature,rtrotter,dim_trott, &
1340 29414 : & tol14,inf_ucvol)
1341 : end do
1342 : end if
1343 :
1344 38 : call timab(609,1,tsec)
1345 29414 : nelectl = sum(rhotry)
1346 38 : call xmpi_sum( nelectl,mpi_enreg%comm_bandfft,ierr)
1347 38 : res_nelectl = inf_ucvol*nelectl - nelect
1348 :
1349 38 : if (res_nelectl /= zero) then
1350 : ! initialisation of fermih
1351 : ! excess of electrons -> smaller fermi
1352 38 : res_nelecth = zero
1353 38 : ii = 1
1354 38 : fermieh = fermie - ten*sign(one,res_nelectl)*temperature
1355 76 : do while(ii<6 .and. res_nelecth*res_nelectl>=0)
1356 38 : fermieh = fermieh - ten*sign(one,res_nelectl)*temperature
1357 38 : call timab(609,2,tsec)
1358 :
1359 38 : if(gputopo) then
1360 : #ifdef HAVE_GPU_CUDA
1361 : call timab(617,1,tsec)
1362 : call density_cuda(npitch,long_tranche,nb_rec,dim_trott,&
1363 : & real(fermieh,cudap),real(temperature,cudap),&
1364 : & real(rtrotter,cudap),real(inf_ucvol,cudap),&
1365 : & real(tol14,cudap),&
1366 : & rhocu)
1367 : rhotry = real(rhocu,dp)
1368 : call timab(617,2,tsec)
1369 : #endif
1370 : else
1371 29414 : do ipointlocal = 1,long_tranche
1372 : call density_rec(a(:,ipointlocal), &
1373 : & b2(:,ipointlocal), &
1374 : & rhotry(ipointlocal), &
1375 : & nb_rec,fermieh,temperature,rtrotter,dim_trott, &
1376 29414 : & tol14,inf_ucvol)
1377 : end do
1378 : end if
1379 38 : call timab(609,1,tsec)
1380 29414 : nelecth = sum(rhotry)
1381 38 : call xmpi_sum( nelecth,mpi_enreg%comm_bandfft ,ierr);
1382 38 : res_nelecth = inf_ucvol*nelecth - nelect
1383 :
1384 38 : if(debug_rec) then
1385 0 : write (msg,'(a,es11.4e2,a,es11.4e2)') ' Fermi energy interval',fermieh,' ',fermiel
1386 0 : call wrtout(std_out,msg,'COLL')
1387 : end if
1388 76 : ii = ii +1
1389 : end do
1390 :
1391 38 : if (res_nelecth*res_nelectl>0) then
1392 0 : write (msg,'(4a)')' fermisolverec : ERROR- ',ch10,&
1393 0 : & ' initial guess for fermi energy doesnt permit to find solutions in solver',ch10
1394 0 : ABI_ERROR(msg)
1395 : end if
1396 :
1397 : ! MAIN LOOP ------------------------------------------------------
1398 143 : main : do nn=1,max_it
1399 : ! fermiem computation
1400 143 : fermiem = 0.5d0*(fermiel+fermieh)
1401 :
1402 : ! nelectm = zero
1403 143 : call timab(609,2,tsec)
1404 :
1405 143 : if(gputopo) then
1406 : #ifdef HAVE_GPU_CUDA
1407 : call timab(617,1,tsec)
1408 : call density_cuda(npitch,long_tranche,nb_rec,dim_trott,&
1409 : & real(fermiem,cudap),real(temperature,cudap),&
1410 : & real(rtrotter,cudap),real(inf_ucvol,cudap),&
1411 : & real(tol14,cudap),&
1412 : & rhocu)
1413 : rhotry = real(rhocu,dp)
1414 : call timab(617,2,tsec)
1415 : #endif
1416 : else
1417 112463 : do ipointlocal = 1,long_tranche
1418 : call density_rec(a(:,ipointlocal), &
1419 : & b2(:,ipointlocal), &
1420 : & rhotry(ipointlocal), &
1421 : & nb_rec,fermiem,temperature,rtrotter,dim_trott, &
1422 112463 : & tol14,inf_ucvol)
1423 : end do
1424 : end if
1425 :
1426 143 : call timab(609,1,tsec)
1427 112463 : nelectm = sum(rhotry)
1428 143 : call xmpi_sum( nelectm,mpi_enreg%comm_bandfft,ierr)
1429 143 : res_nelectm = inf_ucvol*nelectm - nelect
1430 :
1431 : ! new guess
1432 143 : ss = sqrt(res_nelectm**two-res_nelectl*res_nelecth)
1433 143 : fermienew = fermiem + (fermiem-fermiel)*sign(one, res_nelectl-res_nelecth)*res_nelectm/ss
1434 :
1435 143 : call timab(609,2,tsec)
1436 143 : if(gputopo) then
1437 : #ifdef HAVE_GPU_CUDA
1438 : call timab(617,1,tsec)
1439 : call density_cuda(npitch,long_tranche,nb_rec,dim_trott,&
1440 : & real(fermienew,cudap),real(temperature,cudap),&
1441 : & real(rtrotter,cudap),real(inf_ucvol,cudap),&
1442 : & real(tol14,cudap),&
1443 : & rhocu)
1444 : rhotry = real(rhocu,dp)
1445 : call timab(617,2,tsec)
1446 : #endif
1447 : else
1448 112463 : do ipointlocal = 1,long_tranche
1449 : call density_rec(a(:,ipointlocal), &
1450 : & b2(:,ipointlocal), &
1451 : & rhotry(ipointlocal), &
1452 : & nb_rec,fermienew,temperature,rtrotter,dim_trott, &
1453 112463 : & tol14,inf_ucvol)
1454 : end do
1455 : end if
1456 :
1457 143 : call timab(609,1,tsec)
1458 112463 : nelectnew = sum(rhotry)
1459 143 : call xmpi_sum( nelectnew,mpi_enreg%comm_bandfft ,ierr);
1460 143 : res_nelectnew = inf_ucvol*nelectnew - nelect
1461 :
1462 : ! fermiel et fermieh for new iteration
1463 143 : if (sign(res_nelectm,res_nelectnew) /= res_nelectm) then
1464 30 : fermiel = fermiem
1465 30 : res_nelectl = res_nelectm
1466 30 : fermieh = fermienew
1467 30 : res_nelecth = res_nelectnew
1468 113 : else if (sign(res_nelectl,res_nelectnew) /= res_nelectl) then
1469 74 : fermieh = fermienew
1470 74 : res_nelecth = res_nelectnew
1471 39 : else if (sign(res_nelecth,res_nelectnew) /= res_nelecth) then
1472 39 : fermiel = fermienew
1473 39 : res_nelectl = res_nelectnew
1474 : end if
1475 :
1476 : ! are we within the tolerance ?
1477 429 : if ((abs(res_nelectnew) < fermitol).or.(nn == max_it)) then
1478 38 : fermie = fermienew
1479 29414 : rho = rhotry
1480 38 : if(debug_rec) then
1481 0 : write (msg,'(a,es11.4e2,a,i4)')' err, num_iter ', res_nelectnew, ' ',nn
1482 0 : call wrtout(std_out,msg,'COLL')
1483 0 : write(msg,'(a,50a)')' ',('-',ii=1,50)
1484 0 : call wrtout(std_out,msg,'COLL')
1485 : end if
1486 : exit main
1487 : end if
1488 :
1489 : end do main
1490 :
1491 : end if
1492 :
1493 : #ifdef HAVE_GPU_CUDA
1494 : !deallocate array on GPU
1495 : if(gputopo) then
1496 : call dealloc_dens_cuda()
1497 : end if
1498 : call timab(613+swt_tm,1,tsec2) !!--start time-counter: sync gpu-cpu
1499 : call xmpi_barrier(mpi_enreg%comm_bandfft)
1500 : call timab(613+swt_tm,2,tsec2) !!--stop time-counter: sync gpu-cpu
1501 : #endif
1502 :
1503 38 : call timab(609,2,tsec)
1504 38 : end subroutine fermisolverec
1505 : !!***
1506 :
1507 : !!****f* ABINIT/density_rec
1508 : !! NAME
1509 : !! density_rec
1510 : !!
1511 : !! FUNCTION
1512 : !! This routine computes the density using the coefficients corresponding to
1513 : !! continued fraction at a point from a fixed potential.
1514 : !!
1515 : !! INPUTS
1516 : !! coordx, coordy, coordz=coordonnees of the computed point
1517 : !! an, bn2 : coefficient given by density_rec. Input if get_rec_coef=0, output else
1518 : !! nrec=order of density_rec
1519 : !! fermie=fermi energy (Hartree)
1520 : !! tsmear=temperature (Hartree)
1521 : !! rtrotter=real trotter parameter
1522 : !! tol=tolerance criteria for stopping density_rec
1523 : !! inf_ucvol=infinitesimal unit cell volume
1524 : !! dim_trott = max(0,2*trotter-1)
1525 : !!
1526 : !! OUTPUT
1527 : !! rho_out=result of the continued fraction multiplied by a multiplicator
1528 : !!
1529 : !! SIDE EFFECTS
1530 : !!
1531 : !! NOTES
1532 : !! at this time :
1533 : !! - exppot should be replaced by ?
1534 : !! - coord should be replaced by ?
1535 : !! - need a rectangular box (rmet diagonal matrix)
1536 : !!
1537 : !! SOURCE
1538 :
1539 283392 : subroutine density_rec(an,bn2,rho_out,nrec, &
1540 : & fermie,tsmear,rtrotter, &
1541 : & dim_trott,tol,inf_ucvol)
1542 :
1543 : !Arguments -------------------------------
1544 : !scalars
1545 : integer,intent(in) :: nrec
1546 : integer,intent(in) :: dim_trott
1547 : real(dp),intent(in) :: fermie,tol,tsmear,inf_ucvol,rtrotter
1548 : real(dp), intent(out) :: rho_out
1549 : !arrays
1550 : real(dp),intent(in) :: an(0:nrec),bn2(0:nrec)
1551 : !Local variables-------------------------------
1552 : !not used, debugging purpose only
1553 : !for debugging purpose, detailled printing only once for density and ekin
1554 : !scalars
1555 : integer, parameter :: minrec = 3
1556 : integer :: irec
1557 : real(dp) :: beta,mult,prod_b2,error,errold
1558 : real(dp) :: pi_on_rtrotter,twortrotter,exp1,exp2
1559 : complex(dp) :: cinv2rtrotter,coeef_mu,facrec0
1560 : ! character(len=500) :: msg
1561 : !arrays
1562 : real(dp) :: tsec(2)
1563 566784 : complex(dp) :: acc_rho(0:nrec)
1564 566784 : complex(dp) :: D(0:dim_trott),Dold(0:dim_trott)
1565 566784 : complex(dp) :: N(0:dim_trott),Nold(0:dim_trott)
1566 : !**************************************************************************
1567 :
1568 283392 : call timab(605,1,tsec)
1569 :
1570 : !##############################################################
1571 : !--Initialisation of metrics
1572 283392 : mult = two/inf_ucvol !non-spined system
1573 283392 : beta = one/tsmear
1574 :
1575 : !--Variables for optimisation
1576 283392 : pi_on_rtrotter = pi/rtrotter
1577 283392 : twortrotter = two*rtrotter
1578 283392 : exp1 = exp((beta*fermie)/(rtrotter))
1579 283392 : exp2 = exp(beta*fermie/(twortrotter))
1580 283392 : cinv2rtrotter = cmplx(one/twortrotter,zero,dp)
1581 283392 : coeef_mu = cmplx(one/exp2,zero,dp)
1582 :
1583 11619072 : N = czero; D = cone
1584 283392 : facrec0 = cone
1585 11619072 : Nold = czero; Dold = czero
1586 : !--Initialisation of accumulated density
1587 5460480 : acc_rho = czero
1588 : !--Initialisation of estimated error
1589 283392 : prod_b2 = twortrotter/exp1
1590 283392 : errold = zero
1591 :
1592 :
1593 : !##############################################################
1594 : !--Main loop
1595 1353567 : maindo : do irec = 0, nrec
1596 :
1597 : ! ######################################################
1598 : ! --Density computation
1599 : ! !--using the property that: sum_i(bi*c)^2|(z-ai*c)=1/c*sum_i(bi)^2|(z/c-ai)
1600 : ! !and for c =exp(-beta*fermie/(two*rtrotter)
1601 :
1602 : call trottersum(dim_trott,error,prod_b2,pi_on_rtrotter,&
1603 : & facrec0,coeef_mu,exp1,&
1604 : & an(irec),bn2(irec),&
1605 1353567 : & N,D,Nold,Dold)
1606 :
1607 1353567 : if(irec/=nrec .and. irec>=minrec)then
1608 503391 : if((bn2(irec+1)<tol14).or.(mult*error<tol.and.errold<tol)) exit maindo
1609 : end if
1610 1353567 : errold = mult*error
1611 : end do maindo
1612 : !--Accumulated density
1613 5951232 : rho_out = mult*real(cone-sum(N/D)*cinv2rtrotter,dp)
1614 :
1615 283392 : call timab(605,2,tsec)
1616 :
1617 283392 : end subroutine density_rec
1618 : !!***
1619 :
1620 : !!****f* ABINIT/gran_potrec
1621 : !! NAME
1622 : !! gran_potrec
1623 : !!
1624 : !! FUNCTION
1625 : !! This routine computes the local part of the grand-potential at a point using a path integral,
1626 : !! in the recursion method.
1627 : !!
1628 : !! INPUTS
1629 : !! an, bn2 : coefficient given by the recursion.
1630 : !! nrec=order of recursion
1631 : !! trotter=trotter parameter
1632 : !! mult=a multiplicator for computing grand-potential (2 for non-spin-polarized system)
1633 : !! debug_rec=debugging variable
1634 : !! n_pt_integ=points for computation of path integral
1635 : !! xmax= maximum point on the x-axis for integration
1636 : !!
1637 : !! OUTPUT
1638 : !! ene_out=grand-potential at the point
1639 : !! if debug_rec=T then ene_out1,ene_out2,ene_out3,ene_out4 are
1640 : !! the different path branch contriubutions to the grand-potential.
1641 : !! In reality it is not the gren potential but the
1642 : !! grand-potential (omega=-PV) divided by -T
1643 : !!
1644 : !! NOTES
1645 : !! in reality it is not the gren potential but the grand-potential (omega=-PV) divided by -T
1646 : !! at this time :
1647 : !! - mult should be not used
1648 : !! - the routine should be integraly rewrited and use the routine recursion.
1649 : !! - only modified for p /= 0
1650 : !!
1651 : !! SOURCE
1652 :
1653 5238 : subroutine gran_potrec(an,bn2,nrec,trotter,ene_out, mult, &
1654 : & debug_rec,n_pt_integ,xmax,&
1655 : & ene_out1,ene_out2,ene_out3,ene_out4)
1656 :
1657 : !Arguments -------------------------------
1658 : !scalars
1659 : integer,intent(in) :: n_pt_integ,nrec,trotter
1660 : logical,intent(in) :: debug_rec
1661 : real(dp), intent(in) :: mult,xmax
1662 : real(dp),intent(inout) :: ene_out,ene_out1,ene_out2,ene_out3,ene_out4 !vz_i
1663 : !arrays
1664 : real(dp), intent(in) :: an(0:nrec),bn2(0:nrec)
1665 :
1666 : !Local variables-------------------------------
1667 : !scalars
1668 : integer, parameter :: level = 7
1669 : integer, save :: first = 1
1670 : integer :: ii,kk,n_pt_integ_path2
1671 : real(dp) :: epsilon,step,twotrotter,xmin,dr_step
1672 : complex(dp) :: D,Dnew,Dold,N,Nnew,Nold,dz_path,ene_acc,ene_acc1,ene_acc2
1673 : complex(dp) :: ene_acc3,ene_acc4
1674 : complex(dp) :: z_path,delta_calc
1675 : character(len=500) :: message
1676 : !arrays
1677 : real(dp) :: tsec(2)
1678 : ! *************************************************************************
1679 :
1680 :
1681 5238 : call timab(611,1,tsec)
1682 :
1683 : !structured debugging if debug_rec=T : print detailled result the first time we enter gran_potrec
1684 5238 : if(debug_rec .and. first==1)then
1685 0 : write(message,'(a)')' '
1686 0 : call wrtout(std_out,message,'PERS')
1687 0 : write(message,'(a)')' gran_potrec : enter '
1688 0 : call wrtout(std_out,message,'PERS')
1689 0 : write(message,'(a,i8)')'n_pt_integ ' , n_pt_integ
1690 0 : call wrtout(std_out,message,'COLL')
1691 0 : first=0
1692 : end if
1693 :
1694 5238 : ene_out = zero
1695 5238 : ene_acc = czero
1696 5238 : ene_acc1 = czero
1697 5238 : ene_acc2 = czero
1698 5238 : ene_acc3 = czero
1699 5238 : ene_acc4 = czero
1700 :
1701 :
1702 : !path parameters
1703 : !n_pt_integ = 2500
1704 5238 : xmin = -half
1705 5238 : step = (xmax-xmin)/real(n_pt_integ,dp)
1706 5238 : if(trotter==0)then
1707 0 : twotrotter = one
1708 0 : epsilon = .5d-1
1709 : else
1710 5238 : twotrotter = two*real(trotter,dp)
1711 5238 : epsilon = half*sin( pi/twotrotter)
1712 : end if
1713 :
1714 : !xmin = -abs(xmin)**(1.d0/twotrotter)
1715 :
1716 : !####################################################################
1717 : ![xmax + i*epsilon,xmin + i*epsilon]
1718 5238 : dr_step = one/real(n_pt_integ,dp)
1719 5238 : dz_path = -cmplx((xmax-xmin)*dr_step,zero,dp)
1720 534276 : path1: do ii = 0,n_pt_integ
1721 : ! z_path = cmplx(xmin + real(ii,dp)*(xmax-xmin)*dr_step,epsilon,dp)
1722 529038 : z_path = cmplx(xmin,epsilon,dp) - real(ii,dp)*dz_path
1723 529038 : Nold = czero
1724 529038 : Dold = cone
1725 529038 : N = cone
1726 529038 : D = z_path - cmplx(an(0),zero,dp)
1727 :
1728 2116152 : do kk=1,nrec
1729 2116152 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1730 2116152 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1731 :
1732 2116152 : Nold = N
1733 2116152 : Dold = D
1734 2116152 : N = Nnew
1735 2116152 : D = Dnew
1736 :
1737 2645190 : if(kk/=nrec)then
1738 2116152 : if((bn2(kk+1)<tol14))exit
1739 : end if
1740 :
1741 : end do
1742 :
1743 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1744 529038 : delta_calc = func_rec(z_path,twotrotter)* N/D *dz_path
1745 534276 : if(ii==0.or.ii==n_pt_integ)then
1746 10476 : ene_acc = ene_acc + half*delta_calc
1747 10476 : if(debug_rec) ene_acc1 = ene_acc1 + half*delta_calc
1748 : else
1749 518562 : ene_acc = ene_acc + delta_calc
1750 518562 : if(debug_rec) ene_acc1 = ene_acc1 + delta_calc
1751 : end if
1752 : end do path1
1753 :
1754 : !####################################################################
1755 : ![xmin + i*epsilon,xmin]
1756 5238 : if(epsilon/step>4.d0)then
1757 : n_pt_integ_path2 = int(epsilon/step)+1
1758 : else
1759 : n_pt_integ_path2 = 5
1760 : end if
1761 5238 : n_pt_integ_path2 = n_pt_integ
1762 5238 : dr_step = one/real(n_pt_integ_path2,dp)
1763 5238 : dz_path = -cmplx(zero,epsilon*dr_step,dp)
1764 534276 : path2: do ii = 0,n_pt_integ_path2
1765 : ! z_path = cmplx(xmin,real(ii,dp)*epsilon*dr_step,dp)
1766 529038 : z_path = cmplx(xmin,zero,dp)-dz_path*real(ii,dp)
1767 529038 : Nold = czero
1768 529038 : Dold = cone
1769 529038 : N = cone
1770 529038 : D = z_path - cmplx(an(0),zero,dp)
1771 :
1772 2116152 : do kk=1,nrec
1773 2116152 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1774 2116152 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1775 :
1776 2116152 : Nold = N
1777 2116152 : Dold = D
1778 2116152 : N = Nnew
1779 2116152 : D = Dnew
1780 :
1781 2645190 : if(kk/=nrec)then
1782 2116152 : if((bn2(kk+1)<tol14))exit
1783 : end if
1784 :
1785 : end do
1786 :
1787 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1788 529038 : delta_calc = func_rec(z_path,twotrotter)* N/D *dz_path
1789 534276 : if(ii==0.or.ii==n_pt_integ_path2)then
1790 10476 : ene_acc = ene_acc + half*delta_calc
1791 10476 : if(debug_rec) ene_acc3 = ene_acc3 + half*delta_calc
1792 : else
1793 518562 : ene_acc = ene_acc + delta_calc
1794 518562 : if(debug_rec) ene_acc3 = ene_acc3 + delta_calc
1795 : end if
1796 : end do path2
1797 :
1798 :
1799 :
1800 : !####################################################################
1801 : ![xmin,0]
1802 : if(xmin/=czero)then
1803 5238 : dr_step = one/real(n_pt_integ,dp)
1804 5238 : dz_path = cmplx(xmin*dr_step,zero,dp)
1805 529038 : path3: do ii = 1,n_pt_integ !the integrand is 0 at 0
1806 : ! z_path = cmplx(real(ii,dp)*xmin*dr_step,zero,dp)
1807 523800 : z_path = real(ii,dp)*dz_path
1808 :
1809 523800 : Nold = czero
1810 523800 : Dold = cone
1811 523800 : N = cone
1812 523800 : D = z_path - cmplx(an(0),zero,dp)
1813 :
1814 2095200 : do kk=1,nrec
1815 2095200 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1816 2095200 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1817 :
1818 2095200 : Nold = N
1819 2095200 : Dold = D
1820 2095200 : N = Nnew
1821 2095200 : D = Dnew
1822 :
1823 2619000 : if(kk/=nrec)then
1824 2095200 : if((bn2(kk+1)<tol14))exit
1825 : end if
1826 : end do
1827 :
1828 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1829 523800 : delta_calc = func_rec(z_path,twotrotter) * N/D *dz_path
1830 529038 : if(ii==n_pt_integ)then
1831 5238 : ene_acc = ene_acc +half*delta_calc
1832 5238 : if(debug_rec) ene_acc4 = ene_acc4 + half*delta_calc
1833 : else
1834 518562 : ene_acc = ene_acc + delta_calc
1835 518562 : if(debug_rec) ene_acc4 = ene_acc4 +delta_calc
1836 : end if
1837 : end do path3
1838 : end if
1839 :
1840 : !####################################################################
1841 : ![xmax,xmax+i*epsilon]
1842 534276 : dr_step = one/real(n_pt_integ_path2,dp)
1843 534276 : dz_path = cmplx(zero,epsilon*dr_step,dp)
1844 534276 : path4: do ii = 0,n_pt_integ_path2
1845 : ! z_path = cmplx(xmax,real(ii,dp)*epsilon*dr_step,dp)
1846 529038 : z_path = cmplx(xmax,0,dp)+real(ii,dp)*dz_path
1847 :
1848 529038 : Nold = czero
1849 529038 : Dold = cone
1850 529038 : N = cone
1851 529038 : D = z_path - cmplx(an(0),zero,dp)
1852 :
1853 2116152 : do kk=1,nrec
1854 2116152 : Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold
1855 2116152 : Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold
1856 :
1857 2116152 : Nold = N
1858 2116152 : Dold = D
1859 2116152 : N = Nnew
1860 2116152 : D = Dnew
1861 :
1862 2645190 : if(kk/=nrec)then
1863 2116152 : if((bn2(kk+1)<tol14))exit
1864 : end if
1865 :
1866 : end do
1867 :
1868 : ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz
1869 529038 : delta_calc = func_rec(z_path,twotrotter) * N/D *dz_path
1870 534276 : if(ii==0.or.ii==n_pt_integ_path2)then
1871 10476 : ene_acc = ene_acc + half*delta_calc
1872 10476 : if(debug_rec) ene_acc2 = ene_acc2 + half*delta_calc
1873 : else
1874 518562 : ene_acc = ene_acc + delta_calc
1875 518562 : if(debug_rec) ene_acc2 = ene_acc2 + delta_calc
1876 : end if
1877 : end do path4
1878 :
1879 5238 : ene_out = mult*real(ene_acc*cmplx(zero,-piinv,dp),dp)
1880 5238 : if(debug_rec) then
1881 0 : ene_out1 = mult*real(ene_acc1*cmplx(zero,-piinv,dp),dp)
1882 0 : ene_out2 = mult*real(ene_acc2*cmplx(zero,-piinv,dp),dp)
1883 0 : ene_out3 = mult*real(ene_acc3*cmplx(zero,-piinv,dp),dp)
1884 0 : ene_out4 = mult*real(ene_acc4*cmplx(zero,-piinv,dp),dp)
1885 : end if
1886 :
1887 5238 : call timab(611,2,tsec)
1888 :
1889 : contains
1890 :
1891 : !func_rec(z_path,twotrotter) = log(cone+z_path**twotrotter)
1892 :
1893 2110914 : function func_rec(z,x)
1894 :
1895 : complex(dp) :: func_rec
1896 : complex(dp),intent(in) :: z
1897 : real(dp),intent(in) :: x
1898 :
1899 2110914 : func_rec = log(cone+z**x)
1900 :
1901 2110914 : end function func_rec
1902 :
1903 : end subroutine gran_potrec
1904 : !!***
1905 :
1906 : !!****f* ABINIT/nlenergyrec
1907 : !! NAME
1908 : !! nlenergyrec
1909 : !!
1910 : !! FUNCTION
1911 : !! During recursion, it computes the non-local energy
1912 : !!
1913 : !! INPUTS
1914 : !! rset<recursion_type>=contains all recursion parameters
1915 : !! exppot=exponential of -1/tsmear*vtrial (computed only once in vtorhorec)
1916 : !! tsmear=temperature (Hartree)
1917 : !! trotter=trotter parameter
1918 : !! tol=tolerance criteria for stopping recursion_nl
1919 : !! ngfft=information about FFT(dtset%ngfft a priori different from ngfftrec)
1920 : !! mpi_enreg=information about MPI paralelisation
1921 : !! rset<recursion_type> contains all parameter of recursion
1922 : !! typat(natom)=type of pseudo potential associated to any atom
1923 : !! natom=number of atoms
1924 : !!
1925 : !! OUTPUT
1926 : !! enlx=non-local energy
1927 : !!
1928 : !! SIDE EFFECTS
1929 : !!
1930 : !! NOTES
1931 : !!
1932 : !! SOURCE
1933 :
1934 0 : subroutine nlenergyrec(rset,enlx,exppot,ngfft,natom,typat,tsmear,trotter,tol)
1935 :
1936 : !Arguments ------------------------------------
1937 : !Scalar
1938 : integer , intent(in) :: natom,trotter
1939 : real(dp), intent(in) :: tsmear,tol
1940 : type(recursion_type),intent(in) :: rset
1941 : real(dp), intent(out) :: enlx
1942 : !Arrays
1943 : integer , intent(in) :: typat(natom),ngfft(18)
1944 : real(dp), intent(in) :: exppot(0:ngfft(1)*ngfft(2)*ngfft(3)-1)
1945 : !Local variables-------------------------------
1946 : integer :: iatom,jatom
1947 : integer :: ii,ipsp,dim_trott
1948 : integer :: ierr,me_count
1949 : integer :: ilmn,jlmn,ilm,jlm,in,jn,il
1950 : character(len=500) :: msg
1951 : logical :: tronc
1952 : real(dp) :: rho_nl,normali,mult
1953 0 : type(mpi_type):: mpi_loc
1954 : !Arrays
1955 0 : integer :: gcart_loc(3,natom)
1956 : integer :: ngfftrec(3),trasl(3)
1957 : real(dp) :: tsec(2)
1958 0 : real(dp) :: un0(0:rset%nfftrec)
1959 0 : real(dp),pointer :: projec(:,:,:,:,:)
1960 0 : real(dp),allocatable :: exppotloc(:)
1961 0 : real(dp) :: proj_arr(0:rset%ngfftrec(1)-1,0:rset%ngfftrec(2)-1,0:rset%ngfftrec(3)-1)
1962 :
1963 : ! *************************************************************************
1964 :
1965 :
1966 0 : call timab(612,1,tsec) !!--start time-counter: nlenergyrec
1967 :
1968 0 : if(rset%debug)then
1969 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' nlenergyrec : enter'
1970 0 : call wrtout(std_out,msg,'PERS')
1971 : end if
1972 :
1973 0 : write(msg,'(a)')' -- nlenergyrec -----------------------------------'
1974 0 : call wrtout(std_out,msg,'COLL')
1975 :
1976 : !--Initialisation variables
1977 0 : enlx = zero
1978 0 : mult = two !--is twice for non-spinned systems
1979 0 : ngfftrec = rset%ngfftrec(:3)
1980 0 : gcart_loc = rset%inf%gcart
1981 0 : mpi_loc = rset%mpi
1982 0 : me_count = 0
1983 0 : dim_trott = max(0,2*trotter-1)
1984 : nullify(projec)
1985 0 : ABI_MALLOC(projec,(0:rset%ngfftrec(1)-1,0:rset%ngfftrec(2)-1,0:rset%ngfftrec(3)-1,rset%nl%lmnmax,natom))
1986 0 : projec = zero
1987 :
1988 0 : tronc = rset%tronc !--True if troncation is used
1989 0 : if(tronc) then
1990 0 : ABI_MALLOC(exppotloc,(0:rset%nfftrec-1))
1991 : end if
1992 :
1993 :
1994 : !--LOOP ON ATOMS to create projectors-vector
1995 0 : atomloop1: do iatom = 1, natom
1996 0 : ipsp = typat(iatom)
1997 : ! --Aquisition,reshape,translation,rotation of the projectors vector
1998 0 : do ilmn = 1,rset%nl%lmnmax
1999 0 : in = rset%nl%indlmn(3,ilmn,ipsp)
2000 : ! --Projectors vector in 3-composant vector
2001 0 : projec(:,:,:,ilmn,iatom) = reshape(rset%nl%projec(:,ilmn,ipsp),shape=shape(projec(:,:,:,1,1)))
2002 : ! --Moving the projectors vector on the center of the grid
2003 0 : do ii=1,3
2004 0 : projec(:,:,:,ilmn,iatom) = cshift(projec(:,:,:,ilmn,iatom),shift=ngfftrec(ii)/2-gcart_loc(ii,iatom),dim=ii)
2005 : end do
2006 : end do
2007 :
2008 : end do atomloop1
2009 :
2010 :
2011 : !##################################################################
2012 : !--LOOP ON ATOMS (MAIN LOOP)
2013 0 : atomloop: do iatom = 1, natom
2014 0 : ipsp = typat(iatom)
2015 :
2016 : ! --If troncation is present, the considered atom has to be in the
2017 : ! center of the grid so atoms, potential and projectors have to be translated
2018 0 : if(tronc) then
2019 0 : trasl = -rset%inf%gcart(:,iatom)+ngfftrec/2
2020 : ! --Translation of atoms
2021 0 : do jatom=1,natom
2022 0 : gcart_loc(:,jatom) = rset%inf%gcart(:,jatom)+trasl
2023 0 : gcart_loc(:,jatom) = modulo(gcart_loc(:,jatom),ngfft(:3))
2024 : ! --Translation of non-local projectors
2025 0 : do ilmn = 1,rset%nl%lmnmax
2026 0 : projec(:,:,:,ilmn,jatom) = reshape(rset%nl%projec(:,ilmn,typat(jatom)),shape=shape(projec(:,:,:,1,1)))
2027 0 : do ii=1,3
2028 0 : projec(:,:,:,ilmn,jatom) = eoshift(projec(:,:,:,ilmn,jatom),shift=ngfftrec(ii)/2-gcart_loc(ii,jatom),dim=ii)
2029 : end do
2030 : end do
2031 : end do
2032 :
2033 : ! --Translation of the potential
2034 0 : call reshape_pot(trasl,ngfft(1)*ngfft(2)*ngfft(3),rset%nfftrec,ngfft(:3),ngfftrec,exppot,exppotloc)
2035 : end if
2036 :
2037 : ! --Loop on projectors
2038 0 : projloop: do ilmn = 1,rset%nl%lmnmax
2039 0 : me_count = iatom+ilmn*natom-2 !--counter of the number of iteration
2040 : ! --Only the proc me compute
2041 0 : if(mpi_loc%me==mod(me_count,mpi_loc%nproc)) then
2042 0 : ilm = rset%nl%indlmn(4,ilmn,ipsp)
2043 0 : proj_arr = zero
2044 0 : do jlmn = 1,rset%nl%lmnmax
2045 0 : jlm = rset%nl%indlmn(4,jlmn,ipsp)
2046 0 : if(ilm==jlm) then
2047 0 : in = rset%nl%indlmn(3,ilmn,ipsp)
2048 0 : jn = rset%nl%indlmn(3,jlmn,ipsp)
2049 0 : il = rset%nl%indlmn(1,ilmn,ipsp)+1
2050 0 : proj_arr(:,:,:) = proj_arr(:,:,:) + rset%nl%eivec(jn,in,il,ipsp)*projec(:,:,:,jlmn,iatom)
2051 : ! write(std_out,*)'l,m,lm,n,n',il-1,rset%nl%indlmn(2,ilmn,ipsp),ilm,in,jn
2052 : ! write(std_out,*)'eigevectors',rset%nl%eivec(jn,in,il,ipsp)
2053 :
2054 : end if
2055 : end do
2056 :
2057 0 : un0 = pack(proj_arr(:,:,:),mask=.true.)
2058 0 : normali = sum(un0*un0)*rset%inf%ucvol
2059 0 : un0 = (one/sqrt(normali))*un0
2060 :
2061 0 : if(tronc)then
2062 : call recursion_nl(exppotloc,un0,rho_nl,rset,rset%ngfftrec,&
2063 : & tsmear,trotter,dim_trott,tol,typat,&
2064 0 : & natom,projec)
2065 : else
2066 : call recursion_nl(exppot,un0,rho_nl,rset,rset%ngfftrec,&
2067 : & tsmear,trotter,dim_trott,tol,typat,&
2068 0 : & natom,projec)
2069 : end if
2070 :
2071 0 : enlx = enlx+mult*rho_nl*rset%nl%eival(in,il,ipsp)*normali
2072 : end if
2073 :
2074 : end do projloop
2075 : end do atomloop
2076 :
2077 : !--Sum the contribution to the non-local energy computed by any procs
2078 0 : call xmpi_sum(enlx,mpi_loc%comm_bandfft,ierr)
2079 :
2080 0 : if(associated(projec)) then
2081 0 : ABI_FREE(projec)
2082 : end if
2083 0 : if(tronc) then
2084 0 : ABI_FREE(exppotloc)
2085 : end if
2086 :
2087 0 : if(rset%debug)then
2088 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' nlenergyrec : exit'
2089 0 : call wrtout(std_out,msg,'PERS')
2090 : end if
2091 :
2092 0 : call timab(612,2,tsec) !--stop time-counter: nlenergyrec
2093 :
2094 0 : end subroutine nlenergyrec
2095 : !!***
2096 :
2097 :
2098 : !!****f* ABINIT/first_rec
2099 : !! NAME
2100 : !! first_rec
2101 : !!
2102 : !! FUNCTION
2103 : !! When recursion method is used, in the first step this routine
2104 : !! compute some quantities which are used in the rest of the calculation.
2105 : !!
2106 : !! COPYRIGHT
2107 : !! Copyright (C) 2009-2026 ABINIT group (MMancini)
2108 : !! This file is distributed under the terms of the
2109 : !! GNU General Public License, see ~abinit/COPYING
2110 : !! or http://www.gnu.org/copyleft/gpl.txt .
2111 : !!
2112 : !! INPUTS
2113 : !! dtset <type(dataset_type)>=all input variables for this dataset:
2114 : !! | recgratio =fine/coarse grid ratio
2115 : !! | recptrott =trotter parameter
2116 : !! | tsmear =temperature
2117 : !! | recrcut =tut radius in recursion (range of iteration)
2118 : !! | ngfft(18) =FFT grid used as real (fine) grid in recursion
2119 : !! psps <type(pseudopotential_type)>=variables related to pseudo-potentials
2120 : !!
2121 : !! OUTPUT
2122 : !!
2123 : !! SIDE EFFECTS
2124 : !! rset <type(recursion_type)>=variables related to recursion method
2125 : !! | debug<logical> = T if debugging is used
2126 : !! | inf <type(metricrec_type)>=information concerning the infinitesimal metrics
2127 : !! | ngfftrec(18) =truncated (or not, if not ngfftrec=ngfft)FFT grid used as real grid in recursion.
2128 : !! | nfftrec =product(ngfftrec(1:3))
2129 : !! | tronc<logical> = T if truncation is effectively used
2130 : !! | ZT_p = fourier transform of the green_kernel calculated on the fine grid
2131 : !!
2132 : !!
2133 : !! NOTES
2134 : !!
2135 : !! SOURCE
2136 :
2137 6 : subroutine first_rec(dtset,psps,rset)
2138 :
2139 : !Arguments ------------------------------------
2140 : ! scalars
2141 : type(dataset_type),intent(in) :: dtset
2142 : type(pseudopotential_type),intent(in) :: psps
2143 : type(recursion_type),intent(inout) :: rset
2144 : !Local variables-------------------------------
2145 : !scalars
2146 : integer :: nfftrec,trotter,ii,dim_trott
2147 : real(dp) :: tsmear,beta,rtrotter
2148 : character(len=500) :: msg
2149 : !arrays
2150 : integer :: ngfftrec(18)
2151 : real(dp) :: tsec(2)
2152 : #ifdef HAVE_GPU_CUDA
2153 : integer :: max_rec,ierr,testpts,swt_tm
2154 : real(dp) :: rho,tm_ratio
2155 : real(dp) :: time_cu,time_f
2156 : type(recursion_type) :: rset_test
2157 : type(recparall_type) :: parold
2158 : integer :: trasl(3)
2159 : real(dp) :: tsec2(2),tsec3(2)
2160 : real(dp) :: aloc(0,1),b2loc(0,1)
2161 : real(dp) :: dm_projec(0,0,0,1,1)
2162 : real(dp) :: exppot(0:dtset%nfft-1)
2163 : real(dp),allocatable :: exppotloc(:)
2164 : real(cudap),allocatable :: aloc_cu(:),b2loc_cu(:)
2165 : #endif
2166 :
2167 : ! *************************************************************************
2168 :
2169 6 : call timab(601,1,tsec) !!--Start time-counter: initialisation
2170 :
2171 6 : ABI_WARNING("RECURSION")
2172 6 : if(dtset%recgratio>1) then
2173 1 : write(msg,'(a)')'COARSE GRID IS USED'
2174 1 : call wrtout(std_out,msg,'COLL')
2175 : end if
2176 :
2177 : !--Initialisation
2178 6 : trotter = dtset%recptrott !--Trotter parameter
2179 6 : tsmear = dtset%tsmear !--Temperature
2180 6 : beta = one/tsmear !--Inverse of temperature
2181 :
2182 : !--Rewriting the trotter parameter
2183 6 : dim_trott = max(0,2*trotter-1)
2184 6 : rtrotter = max(half,real(trotter,dp))
2185 :
2186 6 : write (msg,'(2a)')ch10,'==== FIRST CYCLE RECURSION ========================='
2187 6 : call wrtout(std_out,msg,'COLL')
2188 :
2189 :
2190 114 : ngfftrec = rset%ngfftrec
2191 6 : nfftrec = rset%nfftrec
2192 : !------------------------------------------------
2193 : !--TRONCATION OF THE BOX: determines new dimensions
2194 : !--Now in InitRec
2195 : !--------------------------------------------------------
2196 : !--DEFINITION PAW VARIABLES COARSE-FINE GRID TO USE TRANSGRID--INGRID FUNCTIONS
2197 : !--Now these variables are defined into gstate by InitRec
2198 :
2199 : !--------------------------------------------------------
2200 : !--COMPUTATION OF THE FOURIER TRANSFORM OF THE GREEN KERNEL (only once)
2201 6 : write (msg,'(a)')' - green kernel calculation -----------------------'
2202 6 : call wrtout(std_out,msg,'COLL')
2203 18 : ABI_MALLOC(rset%ZT_p,(1:2,0: nfftrec-1))
2204 6 : call timab(601,2,tsec)
2205 6 : call green_kernel(rset%ZT_p,rset%inf%rmet,rset%inf%ucvol,rtrotter/beta,rset%mpi,ngfftrec,nfftrec)
2206 6 : call timab(601,1,tsec)
2207 306 : write(msg,'(a,50a)')' ',('-',ii=1,50)
2208 6 : call wrtout(std_out,msg,'COLL')
2209 : !!--end computation of the fourier transform of the Green kernel
2210 :
2211 : !!-----------------------------------
2212 : !!--ROUTINE FOR THE CALCULATION OF THE NON-LOCAL PSEUDO
2213 : !--Now these variables here by Init_nlpspRec
2214 6 : call Init_nlpspRec(four*tsmear*rtrotter,psps,rset%nl,rset%inf,rset%ngfftrec,rset%debug)
2215 :
2216 : !!-----------------------------------
2217 : !--Load distribution on procs when GPU are present
2218 : #if defined HAVE_GPU_CUDA
2219 :
2220 : !--Test timing only if exists GPU and they are not equal to the cpus
2221 : if(rset%tp == 4) then
2222 : parold = rset%par
2223 : ii = 0
2224 : time_f = zero
2225 : time_cu = zero
2226 : call random_number(exppot) ! exppot = one
2227 :
2228 : if(rset%gpudevice == -1) then
2229 : ! --Test CPUS
2230 : swt_tm = 0
2231 : testpts = min(rset%par%npt, 20)
2232 : call timein(tsec2(1),tsec2(2))
2233 : if(rset%tronc) then
2234 : ABI_MALLOC(exppotloc,(0:nfftrec-1))
2235 : do while(ii< testpts)
2236 : trasl = -(/1,2,3/)+ngfftrec(:3)/2
2237 : call reshape_pot(trasl,dtset%nfft,nfftrec,dtset%ngfft(:3),ngfftrec(:3),&
2238 : & exppot,exppotloc)
2239 : call recursion(exppotloc,0,0,0, &
2240 : & aloc, &
2241 : & b2loc, &
2242 : & rho,&
2243 : & 0, rset%efermi,tsmear,rtrotter,dim_trott, &
2244 : & rset%ZT_p, &
2245 : & dtset%rectolden,dtset%typat, &
2246 : & rset%nl,&
2247 : & rset%mpi,nfftrec,ngfftrec,rset%inf,&
2248 : & 6,dtset%natom,dm_projec,0)
2249 : ii=ii+1
2250 : end do
2251 : ABI_FREE(exppotloc)
2252 : else
2253 : do while(ii< testpts)
2254 : call recursion(exppot,0,0,0, &
2255 : & aloc, &
2256 : & b2loc, &
2257 : & rho,&
2258 : & 0, rset%efermi,tsmear,rtrotter,dim_trott, &
2259 : & rset%ZT_p, &
2260 : & dtset%rectolden,dtset%typat, &
2261 : & rset%nl,&
2262 : & rset%mpi,nfftrec,ngfftrec,rset%inf,&
2263 : & 6,dtset%natom,dm_projec,0)
2264 : ii=ii+1
2265 : end do
2266 : end if
2267 : call timein(tsec3(1),tsec3(2))
2268 : time_f = (tsec3(1)-tsec2(1))/real(testpts,dp)
2269 : time_f = time_f*time_f
2270 : else
2271 : ! --Test GPUS
2272 : swt_tm = 1
2273 : rset_test = rset
2274 : rset_test%GPU%par%npt = max(rset%GPU%nptrec,100)
2275 : rset_test%min_nrec = 0
2276 : call get_pt0_pt1(dtset%ngfft(:3),dtset%recgratio,0,&
2277 : & rset_test%GPU%par%npt,rset_test%GPU%par)
2278 :
2279 :
2280 : ABI_MALLOC(aloc_cu,(rset_test%GPU%par%npt))
2281 : ABI_MALLOC(b2loc_cu,(rset_test%GPU%par%npt))
2282 : call timein(tsec2(1),tsec2(2))
2283 : call cudarec(rset_test, exppot,aloc_cu,b2loc_cu,&
2284 : & beta,trotter,dtset%rectolden,dtset%recgratio,dtset%ngfft,max_rec)
2285 : call timein(tsec3(1),tsec3(2))
2286 : ABI_FREE(aloc_cu)
2287 : ABI_FREE(b2loc_cu)
2288 :
2289 : time_cu = (tsec3(1)-tsec2(1))/real(rset_test%GPU%par%npt,dp)
2290 : time_cu = time_cu*time_cu
2291 : end if
2292 :
2293 :
2294 : ! --Get Total Times
2295 : call xmpi_sum(time_f,rset%mpi%comm_bandfft,ierr)
2296 : call xmpi_sum(time_cu,rset%mpi%comm_bandfft,ierr)
2297 :
2298 : ! --Average Total Times
2299 : time_f = sqrt(time_f/real(rset%mpi%nproc-rset%ngpu,dp))
2300 : time_cu = sqrt(time_cu/real(rset%ngpu,dp))
2301 : tm_ratio = time_f/time_cu
2302 :
2303 :
2304 : write(msg,'(3(a25,f10.5,a))')&
2305 : & ' Time for cpu recursion ',time_f,ch10,&
2306 : & ' Time for gpu recursion ',time_cu,ch10,&
2307 : & ' Time ratio ',tm_ratio,ch10
2308 : call wrtout(std_out,msg,'COLL')
2309 :
2310 :
2311 : ! tm_ratio =1.20d2! 0.d0! 1.21d0
2312 : rset%par = parold
2313 : ! --Compute the work-load distribution on devices (gpu,cpu)
2314 : if(tm_ratio>1.5d0 .and. time_cu>zero)then
2315 : rset%load = 1
2316 : call cpu_distribution(dtset%recgratio,rset,dtset%ngfft(:3),tm_ratio,1)
2317 : else
2318 : rset%gpudevice = -1
2319 : end if
2320 : end if
2321 :
2322 : #endif
2323 :
2324 :
2325 : !------------------------------------------------------------
2326 : !--DETERMINING WHICH POINT WILL COMPUTE THAT PROC
2327 : !--Now these variables are defined into gstate by Init_rec
2328 :
2329 6 : write (msg,'(2a)')ch10,'==== END FIRST CYCLE RECURSION ====================='
2330 6 : call wrtout(std_out,msg,'COLL')
2331 6 : call timab(601,2,tsec) !!--stop time-counter: initialisation
2332 :
2333 6 : end subroutine first_rec
2334 : !!***
2335 :
2336 :
2337 : !!****f* ABINIT/green_kernel
2338 : !! NAME
2339 : !! green_kernel
2340 : !!
2341 : !! FUNCTION
2342 : !! this routine compute the fourrier transform of the Green kernel for the
2343 : !! recursion method
2344 : !!
2345 : !! INPUTS
2346 : !! inf_rmet=define the infinitesimal metric : rprimd*(transpose(rprimd)) divided
2347 : !! by the number of discretisation point
2348 : !! inf_ucvol=volume of infinitesimal cell
2349 : !! mult=variance of the Gaussian (=rtrotter/beta)
2350 : !! mpi_enreg=information about MPI parallelization
2351 : !! ngfft=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2352 : !! nfft=total number of fft grid points
2353 : !! debug_rec=debugging variable
2354 : !!
2355 : !! OUTPUT
2356 : !! ZT_p=fourier transforme of the Green kernel
2357 : !!
2358 : !! NOTES
2359 : !! at this time :
2360 : !! - need a rectangular box
2361 : !!
2362 : !! SOURCE
2363 :
2364 :
2365 6 : subroutine green_kernel(ZT_p,inf_rmet,inf_ucvol,mult,mpi_enreg,ngfft,nfft)
2366 :
2367 : !Arguments -------------------------------
2368 : !scalars
2369 : integer,intent(in) :: nfft
2370 : real(dp),intent(in) :: inf_ucvol,mult
2371 : type(MPI_type),intent(in) :: mpi_enreg
2372 : !arrays
2373 : integer,intent(in) :: ngfft(18)
2374 : real(dp),intent(in) :: inf_rmet(3,3)
2375 : real(dp),intent(out) :: ZT_p(1:2,0:nfft-1)
2376 :
2377 : !Local variables-------------------------------
2378 : !scalars
2379 : integer,parameter :: n_green_max=5
2380 : integer :: ii,isign,jj,kk,n_green,xx,yy,zz
2381 : real(dp) :: acc, norme
2382 : character(len=500) :: msg
2383 : !arrays
2384 : real(dp) :: tsec(2)
2385 6 : real(dp),allocatable :: T_p(:)
2386 :
2387 : ! *************************************************************************
2388 :
2389 6 : call timab(603,1,tsec)
2390 :
2391 6 : norme = (mult/pi)**(onehalf)
2392 :
2393 18 : ABI_MALLOC(T_p,(0:nfft-1))
2394 :
2395 : !n_green should be better chosen for non rectangular cell
2396 6 : do xx=1, n_green_max
2397 6 : n_green = xx
2398 : if(exp(-mult*dsq_green(xx*ngfft(1),0,0,inf_rmet))<tol14 &
2399 : & .and. exp(-mult*dsq_green(0,xx*ngfft(2),0,inf_rmet))<tol14 &
2400 6 : & .and. exp(-mult*dsq_green(0,0,xx*ngfft(3),inf_rmet))<tol14 ) exit
2401 : end do
2402 :
2403 8710 : acc = zero
2404 8710 : T_p = zero
2405 70 : do kk = 0,ngfft(3)-1
2406 806 : do jj = 0,ngfft(2)-1
2407 9504 : do ii = 0,ngfft(1)-1
2408 :
2409 26112 : do xx=-n_green,n_green-1
2410 60928 : do yy=-n_green,n_green-1
2411 121856 : do zz=-n_green,n_green-1
2412 :
2413 : T_p(ii+ngfft(1)*jj+ngfft(1)*ngfft(2)*kk) = T_p(ii+ngfft(1)*jj+ngfft(1)*ngfft(2)*kk)+ &
2414 104448 : & exp(-mult*dsq_green(ii+xx*ngfft(1),jj+yy*ngfft(2),kk+zz*ngfft(3),inf_rmet))
2415 :
2416 : end do
2417 : end do
2418 : end do
2419 :
2420 8704 : T_p(ii+ngfft(1)*jj+ngfft(1)*ngfft(2)*kk) = norme*T_p(ii+ngfft(1)*jj+ngfft(1)*ngfft(2)*kk)
2421 9440 : acc = acc + inf_ucvol* T_p(ii+ngfft(1)*jj+ngfft(1)*ngfft(2)*kk)
2422 :
2423 : end do
2424 : end do
2425 : end do
2426 :
2427 8710 : T_p(:)= (one/acc)*T_p(:)
2428 :
2429 : !if(debug_rec)then
2430 : write(msg,'(a,d12.3,2(2a,i8),2(2a,3d12.3),2a,d16.6)')&
2431 6 : & ' on the boundary ', exp(-mult*dsq_green(ngfft(1),0,0,inf_rmet)),ch10, &
2432 8710 : & ' no zero ', count(T_p>tol14),ch10, &
2433 6 : & ' n_green ', n_green,ch10, &
2434 6 : & ' erreur_n_green ', exp(-mult*dsq_green(n_green*ngfft(1),0,0,inf_rmet)), &
2435 6 : & exp(-mult*dsq_green(0,n_green*ngfft(2),0,inf_rmet)), &
2436 6 : & exp(-mult*dsq_green(0,0,n_green*ngfft(3),inf_rmet)),ch10,&
2437 6 : & ' erreur_troncat ', T_p(ngfft(1)/2), &
2438 6 : & T_p(ngfft(1)*(ngfft(2)/2)), &
2439 6 : & T_P(ngfft(1)*ngfft(2)*(ngfft(3)/2)),ch10, &
2440 12 : & ' erreurT_p ',abs(acc-1.d0)
2441 6 : call wrtout(std_out,msg,'COLL')
2442 : !endif
2443 :
2444 :
2445 6 : isign = -1
2446 6 : call fourdp(1,ZT_p,T_p,isign,mpi_enreg,nfft,1,ngfft,0)
2447 :
2448 6 : ABI_FREE(T_p)
2449 :
2450 26118 : ZT_p(:,:) = real(nfft,dp)*ZT_p
2451 :
2452 :
2453 6 : call timab(603,2,tsec)
2454 :
2455 : contains
2456 :
2457 69674 : function dsq_green(ii,jj,kk,inf_rmet)
2458 :
2459 : real(dp) :: dsq_green
2460 : integer,intent(in) :: ii,jj,kk
2461 : real(dp),intent(in) :: inf_rmet(3,3)
2462 : dsq_green= inf_rmet(1,1)*dble(ii**2)&
2463 : & +inf_rmet(2,2)*dble(jj**2)&
2464 : & +inf_rmet(3,3)*dble(kk**2)&
2465 : & +two*(inf_rmet(1,2)*dble(ii*jj)&
2466 : & +inf_rmet(2,3)*dble(jj*kk)&
2467 69674 : & +inf_rmet(3,1)*dble(kk*ii))
2468 69674 : end function dsq_green
2469 :
2470 : end subroutine green_kernel
2471 : !!***
2472 :
2473 :
2474 : !!****f* ABINIT/recursion
2475 : !! NAME
2476 : !! recursion
2477 : !!
2478 : !! FUNCTION
2479 : !! This routine computes the recursion coefficients and the corresponding
2480 : !! continued fraction to get the density at a point from a fixed potential.
2481 : !!
2482 : !! INPUTS
2483 : !! exppot=exponential of -1/tsmear*vtrial (computed only once in vtorhorec)
2484 : !! coordx, coordy, coordz=coordonnees of the computed point
2485 : !! nrec=order of recursion
2486 : !! fermie=fermi energy (Hartree)
2487 : !! tsmear=temperature (Hartree)
2488 : !! dim_trott=dimension of the partial fraction decomposition
2489 : !! rtrotter=trotter parameter (real)
2490 : !! ZT_p=fourier transform of the Green krenel (computed only once in vtorhorec)
2491 : !! typat(:)=type of psp associated to any atom
2492 : !! tol=tolerance criteria for stopping recursion
2493 : !! debug=debugging variable
2494 : !! mpi_enreg=information about MPI paralelisation
2495 : !! nfft=number of points in FFT grid
2496 : !! ngfft=information about FFT
2497 : !! metrec<type(metricrec_type)>=information concerning the infinitesimal metrics
2498 : !! inf_ucvol=infinitesimal unit cell volume
2499 : !! tim_fourdp=time counter for fourdp
2500 : !! natom=number of atoms
2501 : !! projec(ngfftrec(1),ngfftrec(2),ngfftrec(3),lmnmax,natom) is the vector, on the ngfftrec grid containing
2502 : !! the non-lacal projector $Y_{lm}(r-R_A)f_{lk}(r-R_A)
2503 : !! tim= 0 if the time spent in the routine is not taken into account,1 otherwise. For example
2504 : !! when measuring time for loading balancing, we don't want to add the time spent in this to the
2505 : !! total time calculation
2506 : !!
2507 : !! OUTPUT
2508 : !! rho_out=result of the continued fraction multiplied by a multiplicator
2509 : !! an, bn2 : coefficient given by recursion.
2510 : !!
2511 : !! SIDE EFFECTS
2512 : !!
2513 : !! NOTES
2514 : !! at this time :
2515 : !! - exppot should be replaced by ?
2516 : !! - coord should be replaced by ?
2517 : !! - need a rectangular box (rmet diagonal matrix)
2518 : !!
2519 : !! SOURCE
2520 :
2521 17469 : subroutine recursion(exppot,coordx,coordy,coordz,an,bn2,rho_out, &
2522 : & nrec,fermie,tsmear,rtrotter,dim_trott, &
2523 34938 : & ZT_p, tol,typat, &
2524 : & nlrec,mpi_enreg,&
2525 : & nfft,ngfft,metrec,&
2526 17469 : & tim_fourdp,natom,projec,tim)
2527 :
2528 :
2529 : use m_linalg_interfaces
2530 :
2531 : !Arguments -------------------------------
2532 : !scalars
2533 : integer,intent(in) :: coordx,coordy,coordz,nfft,nrec,tim
2534 : integer,intent(in) :: tim_fourdp,natom,dim_trott
2535 : real(dp),intent(in) :: fermie,tol,tsmear,rtrotter
2536 : real(dp), intent(out) :: rho_out
2537 : type(MPI_type),intent(in) :: mpi_enreg
2538 : type(nlpsprec_type),intent(in) :: nlrec
2539 : type(metricrec_type),intent(in) :: metrec
2540 : !arrays
2541 : integer, intent(in) :: ngfft(18)
2542 : integer, intent(in) :: typat(natom)
2543 : real(dp), intent(in) :: ZT_p(1:2, 0:nfft-1)
2544 : real(dp), intent(in) :: exppot(0:nfft-1)
2545 : real(dp), intent(in) :: projec(0:,0:,0:,1:,1:)
2546 : real(dp), intent(out) :: an(0:nrec),bn2(0:nrec)
2547 : !Local variables-------------------------------
2548 : !not used, debugging purpose only
2549 : !for debugging purpose, detailled printing only once for density and ekin
2550 : !scalars
2551 : integer, parameter :: level = 7, minrec = 3
2552 : integer :: irec,isign,timab_id,ii
2553 : real(dp) :: switchimu,switchu
2554 : real(dp) :: bb,beta,mult,prod_b2,error,errold
2555 : real(dp) :: inf_ucvol,pi_on_rtrotter,twortrotter,exp1,exp2
2556 : complex(dp) :: cinv2rtrotter,coeef_mu,facrec0
2557 : ! character(len=500) :: msg
2558 : !arrays
2559 : real(dp) :: tsec(2)
2560 : real(dp) :: inf_tr(3)
2561 34938 : real(dp) :: Zvtempo(1:2, 0:nfft-1)
2562 34938 : real(dp) :: unold(0:nfft-1),vn(0:nfft-1),un(0:nfft-1)
2563 34938 : complex(dp) :: acc_rho(0:nrec)
2564 34938 : complex(dp) :: D(0:dim_trott),Dold(0:dim_trott)
2565 34938 : complex(dp) :: N(0:dim_trott),Nold(0:dim_trott)
2566 : ! *************************************************************************
2567 :
2568 : !--If count time or not
2569 17469 : timab_id = 616; if(tim/=0) timab_id = 606;
2570 :
2571 17469 : call timab(timab_id,1,tsec)
2572 :
2573 : !##############################################################
2574 : !--Initialisation of metrics
2575 17469 : inf_ucvol = metrec%ucvol
2576 : inf_tr = metrec%tr
2577 17469 : mult = two/inf_ucvol !non-spined system
2578 :
2579 17469 : beta = one/tsmear
2580 : !--Variables for optimisation
2581 17469 : pi_on_rtrotter = pi/rtrotter
2582 17469 : twortrotter = two*rtrotter
2583 17469 : exp1 = exp((beta*fermie)/(rtrotter))
2584 17469 : exp2 = exp(beta*fermie/(twortrotter))
2585 17469 : cinv2rtrotter = cmplx(one/twortrotter,zero,dp)
2586 17469 : coeef_mu = cmplx(one/exp2,zero,dp)
2587 :
2588 : !--Initialisation of an,bn,un....
2589 716229 : N = czero; D = cone
2590 17469 : facrec0 = cone
2591 716229 : Nold = czero; Dold = czero
2592 :
2593 568539 : an = zero; bn2 = zero; bn2(0) = one
2594 59761341 : bb = zero; vn = zero; unold = zero
2595 : !--u0 is a Dirac function
2596 29889405 : un = zero
2597 17469 : un(coordx+ngfft(1)*(coordy+ngfft(2)*coordz)) = one/sqrt(inf_ucvol)
2598 :
2599 : !--Initialisation of accumulated density
2600 293004 : acc_rho = czero
2601 : !--Initialisation of estimated error
2602 17469 : prod_b2 = twortrotter/exp1
2603 17469 : errold = zero
2604 :
2605 : !##############################################################
2606 : !--Main loop
2607 69876 : maindo : do irec = 0, nrec
2608 :
2609 : ! --Get an and bn2 coef by the lanczos method
2610 :
2611 : ! --Computation of exp(-beta*V/8*p)*un or exp(-beta*V/4*p)*un
2612 : ! depending on if nl part has to be calculated or not.
2613 119557620 : vn = exppot * un
2614 :
2615 : ! --First Non-local psp contribution: (Id+sum_atom int dr1(E(r,r1))vn(r1))
2616 : ! --Computation of exp(-beta*V_NL/4*p)*vn
2617 69876 : if(nlrec%nlpsp) then
2618 0 : call timab(timab_id,2,tsec)
2619 0 : call vn_nl_rec(vn,natom,typat,ngfft(:3),inf_ucvol,nlrec,projec)
2620 0 : call timab(timab_id,1,tsec)
2621 :
2622 : ! --Computation of exp(-beta*V/8*p)*vn in nonlocal case
2623 0 : vn = exppot * vn
2624 : end if !--End if on nlrec%nlpsp
2625 :
2626 : ! --Convolution with the Green kernel
2627 : ! --FFT of vn
2628 69876 : isign = -1
2629 69876 : call fourdp(1,Zvtempo,vn,isign,mpi_enreg,nfft,1,ngfft,tim_fourdp)
2630 :
2631 : ! --F(T)F(vn)
2632 119557620 : do ii = 0,nfft-1
2633 119487744 : switchu = Zvtempo(1,ii)
2634 119487744 : switchimu = Zvtempo(2,ii)
2635 119487744 : Zvtempo(1,ii) = switchu*ZT_p(1,ii) - switchimu*ZT_p(2,ii)
2636 119557620 : Zvtempo(2,ii) = switchu*ZT_p(2,ii) + switchimu*ZT_p(1,ii)
2637 : end do
2638 :
2639 : ! --F^-1(F(T)F(vn))
2640 69876 : isign = 1
2641 69876 : call fourdp(1,Zvtempo,vn,isign,mpi_enreg,nfft,1,ngfft,tim_fourdp)
2642 :
2643 : ! --Computation of exp(-beta*V/8*p)*un or exp(-beta*V/4*p)*un
2644 : ! depending on if nl part has to be calculated or not.
2645 :
2646 119557620 : vn = inf_ucvol * exppot * vn
2647 :
2648 : ! --Second Non-local psp contribution: (Id+sum_atom E(r,r1))vn
2649 69876 : if(nlrec%nlpsp) then
2650 0 : call timab(timab_id,2,tsec)
2651 0 : call vn_nl_rec(vn,natom,typat,ngfft(:3),inf_ucvol,nlrec,projec)
2652 0 : call timab(timab_id,1,tsec)
2653 :
2654 : ! --Computation of exp(-beta*V/8*p)*vn in nonlocal case
2655 0 : vn = exppot * vn
2656 : end if !--End if on nlrec%nlpsp
2657 :
2658 : ! --Multiplication of a and b2 coef by exp(beta*fermie/(two*rtrotter)) must be done in the continued fraction computation
2659 : ! --Computation of a and b2
2660 69876 : an(irec) = inf_ucvol*ddot(nfft,vn,1,un,1)
2661 :
2662 : ! --an must be positive real
2663 : ! --We must compute bn2 and prepare for the next iteration
2664 69876 : if(irec<nrec)then
2665 119557620 : do ii = 0,nfft-1
2666 119487744 : switchu = un(ii)
2667 119487744 : un(ii) = vn(ii)-an(irec)*un(ii)-bb*unold(ii)
2668 119487744 : unold(ii) = switchu
2669 119557620 : bn2(irec+1) = bn2(irec+1)+inf_ucvol*un(ii)*un(ii)
2670 : end do
2671 69876 : bb = sqrt(bn2(irec+1))
2672 119557620 : un = (one/bb)*un
2673 : end if
2674 :
2675 : ! ######################################################
2676 : ! --Density computation
2677 : ! density computation is done inside the main looping, juste after the calculus of a and b2, in order to make
2678 : ! it possible to stop the recursion at the needed accuracy, without doing more recursion loop than needed -
2679 : ! further developpement
2680 :
2681 : ! !--using the property that: sum_i(bi*c)^2|(z-ai*c)=1/c*sum_i(bi)^2|(z/c-ai)
2682 : ! !and for c =exp(-beta*fermie/(two*rtrotter)
2683 :
2684 :
2685 : call trottersum(dim_trott,error,prod_b2,pi_on_rtrotter,&
2686 : & facrec0,coeef_mu,exp1,&
2687 : & an(irec),bn2(irec),&
2688 69876 : & N,D,Nold,Dold)
2689 :
2690 :
2691 69876 : if(irec/=nrec .and. irec>=minrec)then
2692 17469 : if((bn2(irec+1)<tol14).or.(mult*error<tol.and.errold<tol)) exit
2693 : end if
2694 69876 : errold = mult*error
2695 : end do maindo
2696 : !--Accumulated density
2697 366849 : rho_out = mult*real(cone-sum(N/D)*cinv2rtrotter,dp)
2698 :
2699 :
2700 17469 : call timab(timab_id,2,tsec)
2701 :
2702 17469 : end subroutine recursion
2703 : !!***
2704 :
2705 :
2706 : !!****f* ABINIT/recursion_nl
2707 : !! NAME
2708 : !! recursion_nl
2709 : !!
2710 : !! FUNCTION
2711 : !! Given a $|un>$ vector on the real-space grid this routine calculates
2712 : !! the density in $|un>$ by recursion method.
2713 : !!
2714 : !! INPUTS
2715 : !! exppot=exponential of -1/tsmear*vtrial (computed only once in vtorhorec)
2716 : !! trotter=trotter parameter
2717 : !! dim_trott=dimension of the partial fraction decomposition
2718 : !! tsmear=temperature (Hartree)
2719 : !! tol=tolerance criteria for stopping recursion_nl
2720 : !! ngfft=information about FFT(dtset%ngfft a priori different from ngfftrec)
2721 : !! rset<recursion_type> contains all parameter of recursion
2722 : !! typat(natom)=type of pseudo potential associated to any atom
2723 : !! natom=number of atoms
2724 : !! projec(ngfftrec(1),ngfftrec(2),ngfftrec(3),lmnmax,natom) is the vector, on the ngfftrec grid containing
2725 : !! the non-lacal projector $Y_{lm}(r-R_A)f_{lk}(r-R_A)
2726 : !!
2727 : !! OUTPUT
2728 : !! rho_out=result of the continued fraction multiplied by a multiplicator
2729 : !!
2730 : !! SIDE EFFECTS
2731 : !! un(:,:,:)=initial vector on the grid. it is changed in output
2732 : !!
2733 : !! NOTES
2734 : !! at this time :
2735 : !! - need a rectangular box (rmet diagonal matrix)
2736 : !!
2737 : !! SOURCE
2738 :
2739 0 : subroutine recursion_nl(exppot,un,rho_out,rset,ngfft, &
2740 0 : & tsmear,trotter,dim_trott,tol,typat,&
2741 : & natom,projec)
2742 :
2743 :
2744 : use m_linalg_interfaces
2745 :
2746 : !Arguments -------------------------------
2747 : !scalars
2748 : integer,intent(in) :: trotter,natom,dim_trott
2749 : real(dp),intent(in) :: tol,tsmear
2750 : type(recursion_type),intent(in) :: rset
2751 : real(dp), intent(out) :: rho_out
2752 : !arrays
2753 : integer,intent(in) :: typat(natom),ngfft(18)
2754 : real(dp),intent(in) :: exppot(0:ngfft(1)*ngfft(2)*ngfft(3)-1)
2755 : real(dp),intent(inout) :: un(0:rset%nfftrec-1)
2756 : real(dp),pointer :: projec(:,:,:,:,:)
2757 : !Local variables-------------------------------
2758 : !scalars
2759 : integer, parameter :: minrec = 3
2760 : integer :: irec,isign,ii
2761 : real(dp) :: bb,beta,mult,prod_b2,rtrotter
2762 : real(dp) :: inf_ucvol,pi_on_rtrotter,twortrotter,exp1
2763 : real(dp) :: exp2,error,errold
2764 : real(dp) :: switchu,switchimu
2765 : complex(dp) :: facrec0,cinv2rtrotter,coeef_mu
2766 : character(len=500) :: msg
2767 : type(mpi_type),pointer:: mpi_loc
2768 : !arrays
2769 : real(dp):: tsec(2)
2770 : real(dp):: inf_tr(3)
2771 0 : real(dp):: an(0:rset%min_nrec),bn2(0:rset%min_nrec)
2772 0 : real(dp):: vn(0:rset%nfftrec-1)
2773 0 : real(dp):: unold(0:rset%nfftrec-1)
2774 0 : real(dp):: Zvtempo(1:2,0:rset%nfftrec-1)
2775 0 : complex(dp) :: acc_rho(0:rset%min_nrec)
2776 0 : complex(dp) :: D(0:dim_trott),Dold(0:dim_trott)
2777 0 : complex(dp) :: N(0:dim_trott),Nold(0:dim_trott)
2778 :
2779 : ! *************************************************************************
2780 :
2781 0 : call timab(608,1,tsec) !--start time-counter: recursion_nl
2782 0 : if(rset%debug)then
2783 0 : msg=' '
2784 0 : call wrtout(std_out,msg,'COLL')
2785 : end if
2786 :
2787 : !##############################################################
2788 0 : beta = one/tsmear
2789 :
2790 : !--Rewriting the trotter parameter
2791 0 : rtrotter = max(half,real(trotter,dp))
2792 :
2793 : !--Initialisation of mpi
2794 0 : mpi_loc => rset%mpi
2795 :
2796 : !--Initialisation of metrics
2797 0 : inf_ucvol = rset%inf%ucvol
2798 0 : inf_tr = rset%inf%tr
2799 : mult = one !--In the case of the calculus of the NL-energy
2800 :
2801 : !--Initialisation of an,bn,un....
2802 0 : N = czero; D = cone
2803 0 : facrec0 = cone
2804 0 : Nold = czero; Dold = czero
2805 :
2806 0 : an = zero; bn2 = zero; bn2(0) = one
2807 0 : bb = zero; vn = zero; unold = zero
2808 :
2809 : !--Variables for optimisation
2810 0 : pi_on_rtrotter = pi/rtrotter
2811 0 : twortrotter = two*rtrotter
2812 0 : exp1 = exp((beta*rset%efermi)/(rtrotter))
2813 0 : exp2 = exp(beta*rset%efermi/(twortrotter))
2814 0 : cinv2rtrotter = cmplx(one/twortrotter,zero,dp)
2815 0 : coeef_mu = cmplx(one/exp2,zero,dp)
2816 :
2817 : !--Initialisation of accumulated density
2818 0 : acc_rho = czero
2819 : !--Initialisation of estimated error
2820 0 : prod_b2 = twortrotter/exp1
2821 0 : errold = zero
2822 :
2823 : !##############################################################
2824 : !--Main loop
2825 0 : maindo : do irec = 0, rset%min_nrec
2826 : ! --Get an and bn2 coef by the lanczos method
2827 :
2828 : ! --Computation of exp(-beta*V/8*p)*un
2829 0 : vn = exppot * un
2830 :
2831 : ! --First Non-local psp contribution: (Id+sum_atom E(r,r1))vn
2832 0 : call timab(608,2,tsec)
2833 0 : call vn_nl_rec(vn,natom,typat,rset%ngfftrec(:3),inf_ucvol,rset%nl,projec)
2834 0 : call timab(608,1,tsec)
2835 :
2836 : ! --Computation of exp(-beta*V/8*p)*un
2837 0 : vn = exppot * vn
2838 :
2839 : ! --Convolution with the Green kernel
2840 : ! --FFT of vn
2841 0 : isign = -1
2842 0 : call fourdp(1,Zvtempo,vn,isign,mpi_loc,rset%nfftrec,1,rset%ngfftrec,6)
2843 :
2844 : ! --F(T)F(vn)
2845 0 : do ii = 0,rset%nfftrec-1
2846 0 : switchu = Zvtempo(1,ii)
2847 0 : switchimu = Zvtempo(2,ii)
2848 0 : Zvtempo(1,ii) = switchu*rset%ZT_p(1,ii) - switchimu*rset%ZT_p(2,ii)
2849 0 : Zvtempo(2,ii) = switchu*rset%ZT_p(2,ii) + switchimu*rset%ZT_p(1,ii)
2850 : end do
2851 :
2852 : ! --F^-1(F(T)F(vn))
2853 0 : isign = 1
2854 0 : call fourdp(1,Zvtempo,vn,isign,mpi_loc,rset%nfftrec,1,rset%ngfftrec,6)
2855 :
2856 : ! --Computation of exp(-beta*V/2*p)*vn
2857 0 : vn = inf_ucvol * exppot * vn
2858 :
2859 : ! --Second Non-local psp contribution: (Id+sum_atom E(r,r1))vn
2860 0 : call timab(608,2,tsec)
2861 0 : call vn_nl_rec(vn,natom,typat,rset%ngfftrec(:3),inf_ucvol,rset%nl,projec)
2862 0 : call timab(608,1,tsec)
2863 :
2864 : ! --Computation of exp(-beta*V/8*p)*vn
2865 0 : vn = exppot * vn
2866 :
2867 :
2868 : ! --Multiplication of a and b2 coef by exp(beta*fermie/(2.d0*rtrotter)) must be done in the continued fraction computation
2869 : ! --Computation of a and b2
2870 0 : an(irec) = inf_ucvol*ddot(rset%nfftrec,vn,1,un,1) !--an must be positive real
2871 :
2872 : ! --We must compute bn2 and prepare for the next iteration
2873 0 : if(irec<rset%min_nrec)then
2874 0 : do ii = 0,rset%nfftrec-1
2875 0 : switchu = un(ii)
2876 0 : un(ii) = vn(ii)-an(irec)*un(ii)-bb*unold(ii)
2877 0 : unold(ii) = switchu
2878 0 : bn2(irec+1) = bn2(irec+1)+inf_ucvol*un(ii)*un(ii)
2879 : end do
2880 0 : bb = sqrt(bn2(irec+1))
2881 0 : un = (one/bb)*un
2882 : end if
2883 :
2884 : ! ######################################################
2885 : ! --Density computation
2886 : ! in order to make it possible to stop the recursion_nl at the
2887 : ! needed accuracy, without doing more recursion_nl loop than needed further developpement
2888 :
2889 : call trottersum(dim_trott,error,&
2890 : & prod_b2,pi_on_rtrotter,&
2891 : & facrec0,coeef_mu,exp1,&
2892 : & an(irec),bn2(irec),&
2893 0 : & N,D,Nold,Dold)
2894 :
2895 :
2896 0 : if(irec/=rset%min_nrec .and. irec>=minrec)then
2897 0 : if((bn2(irec+1)<tol14).or.(mult*error<tol.and.errold<tol)) exit
2898 : end if
2899 0 : errold = mult*error
2900 : end do maindo
2901 : !--Accumulated density
2902 0 : rho_out = mult*real(cone-sum(N/D)*cinv2rtrotter,dp)
2903 :
2904 0 : call timab(608,2,tsec) !--stop time-counter: recursion_nl
2905 :
2906 0 : end subroutine recursion_nl
2907 : !!***
2908 :
2909 :
2910 : !!****f* ABINIT/vn_nl_rec
2911 : !! NAME
2912 : !! vn_nl_rec
2913 : !!
2914 : !! FUNCTION
2915 : !! this routine computes the contribution to the vector vn, during
2916 : !! recursion, due to the non-local psp.
2917 : !!
2918 : !! INPUTS
2919 : !! vn(:,:,:)=the vector on the real-space grid.
2920 : !! inf_ucvol=volume of infinitesimal cell
2921 : !! natom=number of atoms
2922 : !! typat(natom)=the type of psps associated to the atoms
2923 : !! ngfftrec(3)=first 3 components of ngfftrec (truncated box, if different from ngfft) for the real-space grid
2924 : !! nlrec<type(nlpsprec_type)> in recursion_type containing information concerning psp
2925 : !! projec(ngfftrec(1),ngfftrec(2),ngfftrec(3),lmnmax,natom) is the vector, on the ngfftrec grid containing
2926 : !! the non-lacal projector $Y_{lm}(r-R_A)f_{lk}(r-R_A)
2927 : !!
2928 : !! OUTPUT
2929 : !! vn_nl(:,:,:)=the non_local contribution to vn
2930 : !!
2931 : !! NOTES
2932 : !!
2933 : !! SOURCE
2934 :
2935 0 : subroutine vn_nl_rec(vn,natom,typat,ngfftrec,inf_ucvol,nlrec,projec)
2936 :
2937 :
2938 : use m_linalg_interfaces
2939 :
2940 : !Arguments -------------------------------
2941 : !scalars
2942 : integer,intent(in) :: natom
2943 : real(dp),intent(in) :: inf_ucvol
2944 : type(nlpsprec_type),intent(in) :: nlrec
2945 : !arrays
2946 : integer,intent(in) :: ngfftrec(3),typat(natom)
2947 : real(dp),intent(in) :: projec(0:,0:,0:,1:,1:)
2948 : real(dp),intent(inout):: vn(0:ngfftrec(1)*ngfftrec(2)*ngfftrec(3)-1)
2949 : !Local variables-------------------------------
2950 : !scalars
2951 : integer :: iatom,nfftrec
2952 : integer :: jlmn,il,in,jn
2953 : integer :: ipsp,ilmn
2954 : integer :: npsp,lmnmax
2955 : real(dp):: vn_nl_loc
2956 : !arrays
2957 0 : real(dp):: vn_nl(0:ngfftrec(1)-1,0:ngfftrec(2)-1,0:ngfftrec(3)-1)
2958 0 : real(dp):: vtempo(0:ngfftrec(1)-1,0:ngfftrec(2)-1,0:ngfftrec(3)-1)
2959 : real(dp):: tsec(2)
2960 : ! *************************************************************************
2961 :
2962 0 : call timab(615,1,tsec)
2963 : !--Initialisation
2964 :
2965 0 : vn_nl = zero
2966 0 : npsp = nlrec%npsp
2967 0 : lmnmax = nlrec%lmnmax
2968 0 : nfftrec = product(ngfftrec)
2969 0 : vtempo(:,:,:) = reshape(source=vn,shape=ngfftrec(:3))
2970 :
2971 : !--Sum_iatom \int dr1 E(r-r_a,r1-r_a)vn(r1) *infucvol
2972 0 : do iatom=1,natom !--Loop on atoms
2973 0 : ipsp = typat(natom)
2974 :
2975 : ! --If psp(typat(iatom)) is local then cycle
2976 0 : if(all(nlrec%pspinfo(:,ipsp)==0)) cycle
2977 :
2978 :
2979 : ! write(std_out,*)'lmnmax',nlrec%lmnmax,lmnmax
2980 :
2981 0 : do ilmn = 1, lmnmax
2982 0 : do jlmn = 1,lmnmax
2983 0 : if(nlrec%indlmn(4,ilmn,ipsp)==nlrec%indlmn(4,jlmn,ipsp)) then
2984 0 : il = 1+nlrec%indlmn(1,jlmn,ipsp)
2985 0 : in = nlrec%indlmn(3,ilmn,ipsp)
2986 0 : jn = nlrec%indlmn(3,jlmn,ipsp)
2987 0 : vn_nl_loc = ddot(nfftrec,projec(:,:,:,jlmn,iatom),1,vtempo,1)
2988 0 : vn_nl = vn_nl+projec(:,:,:,ilmn,iatom)*vn_nl_loc*nlrec%mat_exp_psp_nl(in,jn,il,ipsp)
2989 : end if
2990 : end do
2991 : end do
2992 : end do !--End loop on atoms
2993 0 : vtempo = vtempo + vn_nl*inf_ucvol
2994 :
2995 0 : vn = reshape(source=vtempo,shape=(/nfftrec/))
2996 :
2997 0 : call timab(615,2,tsec)
2998 :
2999 0 : end subroutine vn_nl_rec
3000 : !!***
3001 :
3002 : end module m_vtorhorec
3003 : !!***
|
