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