Branch data Line data Source code
1 : : !!****m* ABINIT/m_wvl_rho
2 : : !! NAME
3 : : !! m_wvl_rho
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2012-2026 ABINIT group (TRangel, DC)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_wvl_rho
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_splines
27 : : use m_errors
28 : : use defs_wvltypes
29 : : use m_sort
30 : : use m_abi2big
31 : : use m_xmpi
32 : : use m_dtset
33 : :
34 : : use defs_abitypes, only : MPI_type
35 : : use m_geometry, only : xred2xcart, metric
36 : : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free
37 : : use m_pawtab, only : pawtab_type
38 : : use m_pawrhoij, only : pawrhoij_type
39 : : use m_drivexc, only : mkdenpos
40 : :
41 : : implicit none
42 : :
43 : : private
44 : : !!***
45 : :
46 : : public :: wvl_initro
47 : : public :: wvl_mkrho
48 : : public :: wvl_prcref
49 : : !!***
50 : :
51 : : contains
52 : : !!***
53 : :
54 : : !!****f* ABINIT/wvl_initro
55 : : !! NAME
56 : : !! wvl_initro
57 : : !!
58 : : !! FUNCTION
59 : : !! FIXME: add description.
60 : : !!
61 : : !! INPUTS
62 : : !! argin(sizein)=description
63 : : !!
64 : : !! OUTPUT
65 : : !! argout(sizeout)=description
66 : : !!
67 : : !! SIDE EFFECTS
68 : : !!
69 : : !! NOTES
70 : : !!
71 : : !! SOURCE
72 : :
73 : 0 : subroutine wvl_initro(&
74 : 0 : & atindx1,geocode,h,me,&
75 : 0 : & natom,nattyp,nfft,nspden,ntypat,&
76 : : & n1,n1i,n2,n2i,n3,&
77 : 0 : & pawrad,pawtab,psppar,&
78 : 0 : & rhor,rprimd,spinat,wvl_den,xc_denpos,xred,zion)
79 : :
80 : : #if defined HAVE_BIGDFT
81 : : use BigDFT_API, only : ELECTRONIC_DENSITY, ext_buffers, ind_positions
82 : : #endif
83 : :
84 : : !Arguments ------------------------------------
85 : : integer,intent(in) :: me,natom,ntypat,nfft,nspden
86 : : integer,intent(in)::n1,n2,n1i,n2i,n3
87 : : real(dp),intent(in) :: h(3)
88 : : type(pawrad_type),intent(in) :: pawrad(ntypat)
89 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
90 : : real(dp),intent(in) :: spinat(3,natom),zion(ntypat)
91 : : real(dp),intent(inout) :: rhor(nfft,nspden)
92 : : real(dp),intent(in) :: xc_denpos
93 : : character(1),intent(in)::geocode
94 : : type(wvl_denspot_type), intent(inout) :: wvl_den
95 : : !arrays
96 : : integer,intent(in) :: atindx1(natom),nattyp(ntypat)
97 : : real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3)
98 : : real(dp),intent(inout)::xred(3,natom)
99 : :
100 : : !Local variables-------------------------------
101 : : #if defined HAVE_BIGDFT
102 : : integer :: ia1,ia2
103 : : integer :: iat,iatm,iatom,iatom_tot,iex,iey,iez,ii,ind
104 : : integer :: ifft,ispden
105 : : integer :: isx,isy,isz,itypat,i1,i2,i3,iwarn,i3s
106 : : integer :: j1,j2,j3,msz
107 : : integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3
108 : : integer :: ncmax,nfgd,nspden_updn,n3pi,shift
109 : : real(dp) :: cutoff,fact,fact0
110 : : real(dp) :: rloc,rr2,rx,ry,rz
111 : : real(dp) :: rshp,r2shp
112 : : real(dp) :: ucvol,xx,yy,zz
113 : : type(pawrad_type)::vale_mesh
114 : : !arrays
115 : : logical :: perx,pery,perz,gox,goy,goz
116 : : real(dp) :: hh(3) !fine grid spacing for wavelets
117 : : real(dp) :: gmet(3,3),gprimd(3,3),rcart(3),rmet(3,3),xcart(3,natom)
118 : : character(len=500) :: message ! to be uncommented, if needed
119 : : !allocatable arrays
120 : : integer,allocatable :: ifftsph_tmp(:),iindex(:)
121 : : real(dp),allocatable:: raux(:),raux2(:)
122 : : real(dp),allocatable:: rr(:)!,rred(:,:)
123 : : #endif
124 : :
125 : : ! *************************************************************************
126 : :
127 : : DBG_ENTER("COLL")
128 : :
129 : : #if defined HAVE_BIGDFT
130 : :
131 : : !PENDING: PARALLELIZATION OVER ATOMS
132 : :
133 : : write(message,'(a,a)') ch10,&
134 : : & ' wvl_initro: Initialize valence density from atomic data by splines'
135 : : call wrtout(std_out,message,'COLL')
136 : :
137 : : !initialize
138 : : rhor(:,:)=zero
139 : :
140 : : if(nspden==4)then
141 : : write(std_out,*)' initro : might work yet for nspden=4 (not checked)'
142 : : write(std_out,*)'spinat',spinat(1:3,1:natom)
143 : : ! stop
144 : : end if
145 : :
146 : : !Check whether the values of spinat are acceptable
147 : : if(nspden==2)then
148 : : do itypat=1,ntypat
149 : : do iat=1,nattyp(itypat)
150 : : iatm=iatm+1;iatom=atindx1(iatm)
151 : : iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom)
152 : :
153 : : if( sqrt(spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2) &
154 : : & > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then
155 : : write(message, '(a,a,a,a,i4,a,a,3es11.4,a,a,a,es11.4)' ) ch10,&
156 : : & ' initro : WARNING - ',ch10,&
157 : : & ' For atom number ',iatom,ch10,&
158 : : & ' input spinat=',spinat(:,iatom),' is larger, in magnitude,',ch10,&
159 : : & ' than zion(ia)=',zion(itypat)
160 : : call wrtout(std_out,message,'COLL')
161 : : call wrtout(ab_out,message,'COLL')
162 : : end if
163 : : end do
164 : : ia1=ia2+1
165 : : end do
166 : : end if
167 : :
168 : : !Fine grid
169 : : hh(:)=0.5d0*h(:)
170 : :
171 : : !mpi:
172 : : !Obtain n3pi, BigDFT quantity:
173 : : n3pi=wvl_den%denspot%dpbox%n3pi
174 : : i3s=wvl_den%denspot%dpbox%nscatterarr(me,3)+1-wvl_den%denspot%dpbox%nscatterarr(me,4)
175 : : shift=n1i*n2i*wvl_den%denspot%dpbox%nscatterarr(me,4)
176 : :
177 : : !Compute xcart from xred
178 : : call xred2xcart(natom,rprimd,xcart,xred)
179 : :
180 : : !Compute metric tensors and ucvol from rprimd
181 : : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
182 : :
183 : : !Conditions for periodicity in the three directions
184 : : perx=(geocode /= 'F')
185 : : pery=(geocode == 'P')
186 : : perz=(geocode /= 'F')
187 : :
188 : : !Compute values of external buffers
189 : : call ext_buffers(perx,nbl1,nbr1)
190 : : call ext_buffers(pery,nbl2,nbr2)
191 : : call ext_buffers(perz,nbl3,nbr3)
192 : :
193 : : iatm=0
194 : : !Big loop on atom types
195 : : do itypat=1,ntypat
196 : : !
197 : : rloc=psppar(0,0,itypat)
198 : : cutoff=10.d0*rloc
199 : :
200 : : ! Create mesh_core object
201 : : ! since tnvale_mesh_size can be bigger than pawrad%mesh_size,
202 : : msz=pawtab(itypat)%tnvale_mesh_size
203 : : call pawrad_init(vale_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,&
204 : : & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep)
205 : : !
206 : : ! Set radius size:
207 : : rshp=vale_mesh%rmax
208 : : r2shp=1.0000001_dp*rshp**2
209 : :
210 : : ! allocate arrays
211 : : if (n3pi > 0) then
212 : : ! sphere: cycle i1,i2,i3
213 : : ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3)
214 : : ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3)
215 : : ! 1+int(1.1* factors are included just for cautioness
216 : : ! circle: cycle only i1 and i2
217 : : ! ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi))
218 : : ! line:
219 : : ncmax=1+int(1.1_dp*rshp/hh(1)*2.d0)
220 : : else
221 : : ncmax=1
222 : : end if
223 : : !
224 : : ABI_MALLOC(ifftsph_tmp,(ncmax))
225 : : ABI_MALLOC(iindex,(ncmax))
226 : : ABI_MALLOC(rr,(ncmax))
227 : : ABI_MALLOC(raux,(ncmax))
228 : : if(nspden==2) then
229 : : ABI_MALLOC(raux2,(ncmax))
230 : : end if
231 : :
232 : : ! Big loop on atoms
233 : : do iat=1,nattyp(itypat)
234 : : iatm=iatm+1;iatom=atindx1(iatm)
235 : : iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom)
236 : :
237 : : ! Spin
238 : : if(nspden==2) then
239 : : fact0=half/zion(itypat)
240 : : fact=fact0*(zion(itypat)+spinat(3,iatom))
241 : : end if
242 : :
243 : : !
244 : : ! Define a "box" around each atom
245 : : rx=xcart(1,iatom_tot)
246 : : ry=xcart(2,iatom_tot)
247 : : rz=xcart(3,iatom_tot)
248 : : !
249 : : isx=floor((rx-cutoff)/hh(1))
250 : : isy=floor((ry-cutoff)/hh(2))
251 : : isz=floor((rz-cutoff)/hh(3))
252 : :
253 : : iex=ceiling((rx+cutoff)/hh(1))
254 : : iey=ceiling((ry+cutoff)/hh(2))
255 : : iez=ceiling((rz+cutoff)/hh(3))
256 : : !
257 : : do i3=isz,iez
258 : : zz=real(i3,kind=8)*hh(3)-rz
259 : : call ind_positions(perz,i3,n3,j3,goz)
260 : : j3=j3+nbl3+1
261 : : !
262 : : do i2=isy,iey
263 : : yy=real(i2,kind=8)*hh(2)-ry
264 : : call ind_positions(pery,i2,n2,j2,goy)
265 : : !
266 : : ! Initialize counters
267 : : nfgd=0
268 : : ! nfgd_r0=0
269 : : !
270 : : do i1=isx,iex
271 : : xx=real(i1,kind=8)*hh(1)-rx
272 : : call ind_positions(perx,i1,n1,j1,gox)
273 : : rr2=xx**2+yy**2+zz**2
274 : : if (j3 >= i3s .and. j3 <= i3s+n3pi-1 .and. goy .and. gox ) then
275 : : !
276 : : if(rr2<=r2shp) then
277 : : if(rr2>tol5) then
278 : : ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
279 : : nfgd=nfgd+1
280 : : rcart=[xx,yy,zz]
281 : : rr(nfgd)=(rr2)**0.5
282 : : ifftsph_tmp(nfgd)=shift+ind
283 : : ! DEBUG
284 : : ! write(itmp,'(i10,3(f13.7,x))')ind,xx+rx,yy+ry,zz+rz
285 : : ! write(itmp,'(6(f13.7,x))')rcart,rred(:,nfgd)
286 : : ! ENDDEBUG
287 : : ! else
288 : : ! ! We save r=0 vectors
289 : : ! ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
290 : : ! ! We reuse the same variable "ifftshp_tmp",
291 : : ! ! but we start from the higher index
292 : : ! nfgd_r0=nfgd_r0+1
293 : : ! ifftsph_tmp(ncmax-nfgd_r0+1)=shift+ind
294 : : end if !rr2>tol5
295 : : end if !rr2<r2shp
296 : : end if !j3..
297 : : end do !i1
298 : :
299 : : ! All of the following could be done inside or outside the loops (i2,i1,i3)
300 : : ! Outside the loops: the memory consuption increases.
301 : : ! Inside the inner loop: the time of calculation increases.
302 : :
303 : : if(nfgd==0) cycle
304 : :
305 : : ! Evaluate spline fit of 1st der of core charge density
306 : : ! from tcoredens(:,2) and tcoredens(:,4)
307 : : do ii=1,nfgd
308 : : iindex(ii)=ii
309 : : end do
310 : : ! write(600,'(i4,x,9999f14.7)')nfgd, rr(1:nfgd)
311 : : call sort_dp(nfgd,rr(1:nfgd),iindex(1:nfgd),tol16)
312 : : call splint(msz,vale_mesh%rad,&
313 : : & pawtab(itypat)%tvalespl(:,1),pawtab(itypat)%tvalespl(:,2),&
314 : : & nfgd,rr(1:nfgd),raux(1:nfgd))
315 : :
316 : :
317 : : ! Accumulate contributions to valence density on the entire cell
318 : : rhor(ifftsph_tmp(1:nfgd),1)=rhor(ifftsph_tmp(1:nfgd),1)+raux(iindex(1:nfgd))
319 : :
320 : : if(nspden==2) then
321 : : raux2(1:nfgd)=raux(iindex(1:nfgd))*fact
322 : : rhor(ifftsph_tmp(1:nfgd),2)=rhor(ifftsph_tmp(1:nfgd),1)+raux2(1:nfgd)
323 : : end if
324 : : ! DEBUG
325 : : ! do ii=1,msz
326 : : ! write(itmp,'(2(f15.7,1x))')vale_mesh%rad(ii),pawtab(itypat)%tvalespl(ii,1)
327 : : ! end do
328 : : ! do ii=1,nfgd
329 : : ! write(itmp,'(2i10)')ii,iindex(ii)
330 : : ! write(itmp,'(2(f15.7,1x))')rr(iindex(ii)),rhor(ifftsph_tmp(ii),1)!,raux(iindex(ii))
331 : : ! end do
332 : : ! END DEBUG
333 : :
334 : : end do !i2
335 : : end do !i1
336 : : end do !iat
337 : :
338 : : ! Deallocate
339 : : call pawrad_free(vale_mesh)
340 : : ABI_FREE(ifftsph_tmp)
341 : : ABI_FREE(iindex)
342 : : ABI_FREE(rr)
343 : : ABI_FREE(raux)
344 : : if(nspden==2) then
345 : : ABI_FREE(raux2)
346 : : end if
347 : :
348 : : end do !itypat
349 : :
350 : : !nspden_updn: 1 for non-polarized, 2 for polarized
351 : : nspden_updn=min(nspden,2)
352 : :
353 : : !Make the density positive everywhere
354 : : call mkdenpos(iwarn,nfft,nspden_updn,1,rhor(:,1:nspden_updn),xc_denpos)
355 : :
356 : : !There seems to be a bug in the intel11 compiler
357 : : !rhor = reshape(wvl_den%denspot%rhov, shape(rhor))
358 : : do ispden=1,nspden
359 : : do ifft=1,nfft
360 : : ii=ifft+nfft*(ispden-1)
361 : : ! rhor(ifft,ispden)=wvl_den%denspot%rhov(ii)
362 : : wvl_den%denspot%rhov(ii)=rhor(ifft,ispden)
363 : : end do
364 : : end do
365 : : wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY
366 : : write(message, '(a,a,a,a)' ) ch10, ' wvl_initro : but why are you copying me :..o('
367 : : call wrtout(std_out,message,'COLL')
368 : :
369 : :
370 : : #else
371 : 0 : BIGDFT_NOTENABLED_ERROR()
372 : : if (.false.) write(std_out,*) me,natom,ntypat,nfft,nspden,n1,n2,n1i,n2i,n3,h(1),&
373 : : & pawrad(1)%mesh_size,pawtab(1)%mesh_size,spinat(1,1),zion(1),rhor(1,1),xc_denpos,&
374 : : & geocode,wvl_den%symObj,atindx1(1),nattyp(1),psppar(1,1,1),rprimd(1,1),xred(1,1)
375 : : #endif
376 : :
377 : : DBG_EXIT("COLL")
378 : :
379 : 0 : end subroutine wvl_initro
380 : : !!***
381 : :
382 : : !!****f* ABINIT/wvl_mkrho
383 : : !! NAME
384 : : !! wvl_mkrho
385 : : !!
386 : : !! FUNCTION
387 : : !! This method is just a wrapper around the BigDFT routine to compute the
388 : : !! density from the wavefunctions.
389 : : !!
390 : : !! INPUTS
391 : : !! dtset <type(dataset_type)>=input variables.
392 : : !! mpi_enreg=information about MPI parallelization
393 : : !! occ(dtset%mband)=occupation numbers.
394 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
395 : : !! wvl_wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.
396 : : !!
397 : : !! OUTPUT
398 : : !! rhor(dtset%nfft)=electron density in r space
399 : : !!
400 : : !! SIDE EFFECTS
401 : : !! proj <type(wvl_projector_type)>=projectors information for wavelets.
402 : : !! | proj(OUT)=computed projectors.
403 : : !!
404 : : !! SOURCE
405 : :
406 : 0 : subroutine wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den)
407 : :
408 : : #if defined HAVE_BIGDFT
409 : : use BigDFT_API, only : sumrho, symmetry_data, ELECTRONIC_DENSITY, communicate_density
410 : : #endif
411 : :
412 : : !Arguments -------------------------------
413 : : !scalars
414 : : type(MPI_type),intent(in) :: mpi_enreg
415 : : type(dataset_type),intent(in) :: dtset
416 : : type(wvl_wf_type),intent(inout) :: wvl_wfs
417 : : type(wvl_denspot_type), intent(inout) :: wvl_den
418 : : !arrays
419 : : real(dp),intent(inout) :: rhor(dtset%nfft,dtset%nspden)
420 : : integer, target, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2, &
421 : : & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
422 : : real(dp), target, intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym), &
423 : : & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
424 : :
425 : : !Local variables-------------------------------
426 : : #if defined HAVE_BIGDFT
427 : : !scalars
428 : : character(len=500) :: message
429 : : integer :: comm,me,nproc
430 : : type(symmetry_data) :: sym
431 : : !for debugging:
432 : : !integer::ifile,ierr
433 : : #endif
434 : :
435 : : ! *************************************************************************
436 : :
437 : : DBG_ENTER("COLL")
438 : :
439 : : #if defined HAVE_BIGDFT
440 : : comm=mpi_enreg%comm_wvl
441 : : me=xmpi_comm_rank(comm)
442 : : nproc=xmpi_comm_size(comm)
443 : :
444 : : sym%symObj = wvl_den%symObj
445 : : sym%irrzon => irrzon
446 : : sym%phnons => phnons
447 : :
448 : : call sumrho(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs,wvl_wfs%ks%Lzd,&
449 : : & wvl_wfs%GPU,sym,wvl_den%denspot%rhod,wvl_den%denspot%xc,&
450 : : & wvl_wfs%ks%psi,wvl_den%denspot%rho_psi)
451 : :
452 : : call communicate_density(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs%nspin,&
453 : : & wvl_den%denspot%rhod,wvl_den%denspot%rho_psi,wvl_den%denspot%rhov,.false.)
454 : :
455 : : wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY
456 : : write(message, '(a,a,a,a)' ) ch10, ' wvl_mkrho : but why are you copying me :..o('
457 : : call wrtout(std_out,message,'COLL')
458 : :
459 : : call wvl_rho_abi2big(2,rhor,wvl_den)
460 : :
461 : : #else
462 : 0 : BIGDFT_NOTENABLED_ERROR()
463 : : if (.false.) write(std_out,*) mpi_enreg%me,dtset%nstep,wvl_wfs%ks,wvl_den%symObj,&
464 : : & rhor(1,1),irrzon(1,1,1),phnons(1,1,1)
465 : : #endif
466 : :
467 : : DBG_EXIT("COLL")
468 : :
469 : 0 : end subroutine wvl_mkrho
470 : : !!***
471 : :
472 : : !!****f* ABINIT/wvl_prcref
473 : : !! NAME
474 : : !! wvl_prcref
475 : : !!
476 : : !! FUNCTION
477 : : !! FIXME: add description.
478 : : !!
479 : : !! INPUTS
480 : : !! argin(sizein)=description
481 : : !!
482 : : !! OUTPUT
483 : : !! argout(sizeout)=description
484 : : !!
485 : : !! SIDE EFFECTS
486 : : !!
487 : : !! NOTES
488 : : !!
489 : : !! SOURCE
490 : :
491 : 0 : subroutine wvl_prcref(dielar,iprcel,my_natom,nfftprc,npawmix,nspden,pawrhoij,&
492 : 0 : & rhoijrespc,usepaw,vresid,vrespc)
493 : :
494 : : !Arguments ------------------------------------
495 : : integer , intent(in) :: iprcel,nfftprc,my_natom,npawmix,nspden,usepaw
496 : : real(dp), intent(in) :: dielar(7)
497 : : real(dp), intent(in) :: vresid(nfftprc,nspden)
498 : : real(dp),intent(out) :: rhoijrespc(npawmix)
499 : : real(dp),intent(out) :: vrespc(nfftprc,nspden)
500 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
501 : :
502 : : !Local variables-------------------------------
503 : : integer :: iatom,index,ispden,klmn,kmix
504 : : real(dp):: diemix,diemixmag,mixfac_eff
505 : : character(len=500) :: message ! to be uncommented, if needed
506 : :
507 : : ! *************************************************************************
508 : :
509 : : DBG_ENTER("COLL")
510 : :
511 : : !PENDING:
512 : : !optres==1, or 0, for density or potential mixing.
513 : : !for potential mixing, we have to average over spins.
514 : : !check prcref.F90 and moddiel.F90
515 : :
516 : :
517 : : #if defined HAVE_BIGDFT
518 : : #endif
519 : :
520 [ # # ]: 0 : if(iprcel .ne. 0) then
521 : : write(message, '(a,i3,a,a,a,a)' )&
522 : 0 : & ' From the calling routine, iprcel=',iprcel,ch10,&
523 : 0 : & ' For wavelets, the only allowed value is 0.',ch10,&
524 : 0 : & ' Action : correct your input file.'
525 : 0 : ABI_ERROR(message)
526 : : end if
527 : :
528 : : !call wvl_moddiel !PENDING
529 : 0 : diemix=dielar(4)
530 : : !dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ;
531 : 0 : diemixmag=abs(dielar(7))
532 [ # # ]: 0 : vrespc(:,1)=diemix*vresid(:,1)
533 [ # # # # : 0 : if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
# # ]
534 : :
535 : : !3) PAW only : precondition the rhoij quantities (augmentation
536 : : !occupancies) residuals. Use a simple preconditionning
537 : : !with the same mixing factor as the model dielectric function.
538 : :
539 [ # # # # ]: 0 : if (usepaw==1.and.my_natom>0) then
540 [ # # ]: 0 : ABI_CHECK(pawrhoij(1)%qphase==1,'wvl_prcref: not available with qphase=1!')
541 : : ! mixfac=dielar(4);mixfacmag=abs(dielar(7))
542 [ # # ]: 0 : if (pawrhoij(1)%cplex_rhoij==1) then
543 : : index=0
544 [ # # ]: 0 : do iatom=1,my_natom
545 [ # # ]: 0 : do ispden=1,pawrhoij(iatom)%nspden
546 [ # # ]: 0 : mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag
547 [ # # ]: 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
548 : 0 : index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
549 : 0 : rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
550 : : end do
551 : : end do
552 : : end do
553 : : else
554 : : index=-1
555 [ # # ]: 0 : do iatom=1,my_natom
556 [ # # ]: 0 : do ispden=1,pawrhoij(iatom)%nspden
557 [ # # ]: 0 : mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag
558 [ # # ]: 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
559 : 0 : index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
560 [ # # ]: 0 : rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
561 : : end do
562 : : end do
563 : : end do
564 : : end if
565 : : end if
566 : :
567 : :
568 : :
569 : : DBG_EXIT("COLL")
570 : :
571 : 0 : end subroutine wvl_prcref
572 : : !!***
573 : :
574 : : end module m_wvl_rho
575 : : !!***
|
