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