TABLE OF CONTENTS


ABINIT/m_wvl_rho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_wvl_rho

FUNCTION

COPYRIGHT

  Copyright (C) 2012-2018 ABINIT group (TRangel, DC)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_wvl_rho
28 
29  use defs_basis
30  use m_abicore
31  use m_splines
32  use m_errors
33  use defs_wvltypes
34  use m_sort
35  use defs_abitypes
36  use m_abi2big
37  use m_xmpi
38 
39  use m_geometry, only : xred2xcart, metric
40  use m_pawrad,   only : pawrad_type, pawrad_init, pawrad_free
41  use m_pawtab,   only : pawtab_type
42  use m_pawrhoij, only : pawrhoij_type
43  use m_drivexc,  only : mkdenpos
44 
45  implicit none
46 
47  private

ABINIT/wvl_initro [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_initro

FUNCTION

  FIXME: add description.

INPUTS

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

PARENTS

      gstate

CHILDREN

      ext_buffers,ind_positions,metric,mkdenpos,pawrad_free,pawrad_init
      sort_dp,splint,wrtout,xred2xcart

SOURCE

 84 subroutine wvl_initro(&
 85 & atindx1,geocode,h,me,&
 86 & natom,nattyp,nfft,nspden,ntypat,&
 87 & n1,n1i,n2,n2i,n3,&
 88 & pawrad,pawtab,psppar,&
 89 & rhor,rprimd,spinat,wvl_den,xc_denpos,xred,zion)
 90 
 91 #if defined HAVE_BIGDFT
 92   use BigDFT_API, only : ELECTRONIC_DENSITY, ext_buffers, ind_positions
 93 #endif
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'wvl_initro'
 99 !End of the abilint section
100 
101  implicit none
102 
103 !Arguments ------------------------------------
104  integer,intent(in) :: me,natom,ntypat,nfft,nspden
105  integer,intent(in)::n1,n2,n1i,n2i,n3
106  real(dp),intent(in) :: h(3)
107  type(pawrad_type),intent(in) :: pawrad(ntypat)
108  type(pawtab_type),intent(in) :: pawtab(ntypat)
109  real(dp),intent(in) :: spinat(3,natom),zion(ntypat)
110  real(dp),intent(inout) :: rhor(nfft,nspden)
111  real(dp),intent(in) :: xc_denpos
112  character(1),intent(in)::geocode
113  type(wvl_denspot_type), intent(inout) :: wvl_den
114 !arrays
115  integer,intent(in) :: atindx1(natom),nattyp(ntypat)
116  real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3)
117  real(dp),intent(inout)::xred(3,natom)
118 
119 !Local variables-------------------------------
120 #if defined HAVE_BIGDFT
121  integer  :: ia1,ia2
122  integer  :: iat,iatm,iatom,iatom_tot,iex,iey,iez,ii,ind
123  integer  :: ifft,ispden
124  integer  :: isx,isy,isz,itypat,i1,i2,i3,iwarn,i3s
125  integer  :: j1,j2,j3,msz
126  integer  :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3
127  integer  :: ncmax,nfgd,nspden_updn,n3pi,shift
128  real(dp) :: cutoff,fact,fact0
129  real(dp) :: rloc,rr2,rx,ry,rz
130  real(dp) :: rshp,r2shp
131  real(dp) :: ucvol,xx,yy,zz
132  type(pawrad_type)::vale_mesh
133 !arrays
134  logical :: perx,pery,perz,gox,goy,goz
135  real(dp) :: hh(3) !fine grid spacing for wavelets
136  real(dp) :: gmet(3,3),gprimd(3,3),rcart(3),rmet(3,3),xcart(3,natom)
137  character(len=500) :: message                   ! to be uncommented, if needed
138 !allocatable arrays
139  integer,allocatable :: ifftsph_tmp(:),iindex(:)
140  real(dp),allocatable:: raux(:),raux2(:)
141  real(dp),allocatable:: rr(:)!,rred(:,:)
142 #endif
143 
144 ! *************************************************************************
145 
146  DBG_ENTER("COLL")
147 
148 #if defined HAVE_BIGDFT
149 
150 !PENDING: PARALLELIZATION OVER ATOMS
151 
152  write(message,'(a,a)') ch10,&
153 & ' wvl_initro: Initialize valence density from atomic data by splines'
154  call wrtout(std_out,message,'COLL')
155 
156 !initialize
157  rhor(:,:)=zero
158 
159  if(nspden==4)then
160    write(std_out,*)' initro : might work yet for nspden=4 (not checked)'
161    write(std_out,*)'spinat',spinat(1:3,1:natom)
162 !  stop
163  end if
164 
165 !Check whether the values of spinat are acceptable
166  if(nspden==2)then
167    do itypat=1,ntypat
168      do iat=1,nattyp(itypat)
169        iatm=iatm+1;iatom=atindx1(iatm)
170        iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom)
171 
172        if( sqrt(spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2) &
173 &       > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then
174          write(message, '(a,a,a,a,i4,a,a,3es11.4,a,a,a,es11.4)' ) ch10,&
175 &         ' initro : WARNING - ',ch10,&
176 &         '  For atom number ',iatom,ch10,&
177 &         '  input spinat=',spinat(:,iatom),'  is larger, in magnitude,',ch10,&
178 &         '  than zion(ia)=',zion(itypat)
179          call wrtout(std_out,message,'COLL')
180          call wrtout(ab_out,message,'COLL')
181        end if
182      end do
183      ia1=ia2+1
184    end do
185  end if
186 
187 !Fine grid
188  hh(:)=0.5d0*h(:)
189 
190 !mpi:
191 !Obtain n3pi, BigDFT quantity:
192  n3pi=wvl_den%denspot%dpbox%n3pi
193  i3s=wvl_den%denspot%dpbox%nscatterarr(me,3)+1-wvl_den%denspot%dpbox%nscatterarr(me,4)
194  shift=n1i*n2i*wvl_den%denspot%dpbox%nscatterarr(me,4)
195 
196 !Compute xcart from xred
197  call xred2xcart(natom,rprimd,xcart,xred)
198 
199 !Compute metric tensors and ucvol from rprimd
200  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
201 
202 !Conditions for periodicity in the three directions
203  perx=(geocode /= 'F')
204  pery=(geocode == 'P')
205  perz=(geocode /= 'F')
206 
207 !Compute values of external buffers
208  call ext_buffers(perx,nbl1,nbr1)
209  call ext_buffers(pery,nbl2,nbr2)
210  call ext_buffers(perz,nbl3,nbr3)
211 
212  iatm=0
213 !Big loop on atom types
214  do itypat=1,ntypat
215 !
216    rloc=psppar(0,0,itypat)
217    cutoff=10.d0*rloc
218 
219 !  Create mesh_core object
220 !  since tnvale_mesh_size can be bigger than pawrad%mesh_size,
221    msz=pawtab(itypat)%tnvale_mesh_size
222    call pawrad_init(vale_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,&
223 &   rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep)
224 !
225 !  Set radius size:
226    rshp=vale_mesh%rmax
227    r2shp=1.0000001_dp*rshp**2
228 
229 !  allocate arrays
230    if (n3pi > 0) then
231 !    sphere: cycle i1,i2,i3
232 !    ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3)
233 !    ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3)
234 !    1+int(1.1* factors are included just for cautioness
235 !    circle: cycle only i1 and i2
236 !    ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi))
237 !    line:
238      ncmax=1+int(1.1_dp*rshp/hh(1)*2.d0)
239    else
240      ncmax=1
241    end if
242 !
243    ABI_ALLOCATE(ifftsph_tmp,(ncmax))
244    ABI_ALLOCATE(iindex,(ncmax))
245    ABI_ALLOCATE(rr,(ncmax))
246    ABI_ALLOCATE(raux,(ncmax))
247    if(nspden==2) then
248      ABI_ALLOCATE(raux2,(ncmax))
249    end if
250 
251 !  Big loop on atoms
252    do iat=1,nattyp(itypat)
253      iatm=iatm+1;iatom=atindx1(iatm)
254      iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom)
255 
256 !    Spin
257      if(nspden==2) then
258        fact0=half/zion(itypat)
259        fact=fact0*(zion(itypat)+spinat(3,iatom))
260      end if
261 
262 !
263 !    Define a "box" around each atom
264      rx=xcart(1,iatom_tot)
265      ry=xcart(2,iatom_tot)
266      rz=xcart(3,iatom_tot)
267 !
268      isx=floor((rx-cutoff)/hh(1))
269      isy=floor((ry-cutoff)/hh(2))
270      isz=floor((rz-cutoff)/hh(3))
271 
272      iex=ceiling((rx+cutoff)/hh(1))
273      iey=ceiling((ry+cutoff)/hh(2))
274      iez=ceiling((rz+cutoff)/hh(3))
275 !
276      do i3=isz,iez
277        zz=real(i3,kind=8)*hh(3)-rz
278        call ind_positions(perz,i3,n3,j3,goz)
279        j3=j3+nbl3+1
280 !
281        do i2=isy,iey
282          yy=real(i2,kind=8)*hh(2)-ry
283          call ind_positions(pery,i2,n2,j2,goy)
284 !
285 !        Initialize counters
286          nfgd=0
287 !        nfgd_r0=0
288 !
289          do i1=isx,iex
290            xx=real(i1,kind=8)*hh(1)-rx
291            call ind_positions(perx,i1,n1,j1,gox)
292            rr2=xx**2+yy**2+zz**2
293            if (j3 >= i3s .and. j3 <= i3s+n3pi-1  .and. goy  .and. gox ) then
294 !
295              if(rr2<=r2shp) then
296                if(rr2>tol5) then
297                  ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
298                  nfgd=nfgd+1
299                  rcart=[xx,yy,zz]
300                  rr(nfgd)=(rr2)**0.5
301                  ifftsph_tmp(nfgd)=shift+ind
302 !                DEBUG
303 !                write(itmp,'(i10,3(f13.7,x))')ind,xx+rx,yy+ry,zz+rz
304 !                write(itmp,'(6(f13.7,x))')rcart,rred(:,nfgd)
305 !                ENDDEBUG
306 !                else
307 !                !              We save r=0 vectors
308 !                ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
309 !                !              We reuse the same variable "ifftshp_tmp",
310 !                !              but we start from the higher index
311 !                nfgd_r0=nfgd_r0+1
312 !                ifftsph_tmp(ncmax-nfgd_r0+1)=shift+ind
313                end if !rr2>tol5
314              end if !rr2<r2shp
315            end if !j3..
316          end do !i1
317 
318 !        All of the following  could be done inside or outside the loops (i2,i1,i3)
319 !        Outside the loops: the memory consuption increases.
320 !        Inside the inner loop: the time of calculation increases.
321 
322          if(nfgd==0)      cycle
323 
324 !        Evaluate spline fit of 1st der of core charge density
325 !        from tcoredens(:,2) and tcoredens(:,4)
326          do ii=1,nfgd
327            iindex(ii)=ii
328          end do
329 !        write(600,'(i4,x,9999f14.7)')nfgd, rr(1:nfgd)
330          call sort_dp(nfgd,rr(1:nfgd),iindex(1:nfgd),tol16)
331          call splint(msz,vale_mesh%rad,&
332 &         pawtab(itypat)%tvalespl(:,1),pawtab(itypat)%tvalespl(:,2),&
333 &         nfgd,rr(1:nfgd),raux(1:nfgd))
334 
335 
336 !        Accumulate contributions to valence density on the entire cell
337          rhor(ifftsph_tmp(1:nfgd),1)=rhor(ifftsph_tmp(1:nfgd),1)+raux(iindex(1:nfgd))
338 
339          if(nspden==2) then
340            raux2(1:nfgd)=raux(iindex(1:nfgd))*fact
341            rhor(ifftsph_tmp(1:nfgd),2)=rhor(ifftsph_tmp(1:nfgd),1)+raux2(1:nfgd)
342          end if
343 !        DEBUG
344 !        do ii=1,msz
345 !        write(itmp,'(2(f15.7,1x))')vale_mesh%rad(ii),pawtab(itypat)%tvalespl(ii,1)
346 !        end do
347 !        do ii=1,nfgd
348 !        write(itmp,'(2i10)')ii,iindex(ii)
349 !        write(itmp,'(2(f15.7,1x))')rr(iindex(ii)),rhor(ifftsph_tmp(ii),1)!,raux(iindex(ii))
350 !        end do
351 !        END DEBUG
352 
353        end do !i2
354      end do !i1
355    end do !iat
356 
357 !  Deallocate
358    call pawrad_free(vale_mesh)
359    ABI_DEALLOCATE(ifftsph_tmp)
360    ABI_DEALLOCATE(iindex)
361    ABI_DEALLOCATE(rr)
362    ABI_DEALLOCATE(raux)
363    if(nspden==2) then
364      ABI_DEALLOCATE(raux2)
365    end if
366 
367  end do !itypat
368 
369 !nspden_updn: 1 for non-polarized, 2 for polarized
370  nspden_updn=min(nspden,2)
371 
372 !Make the density positive everywhere
373  call mkdenpos(iwarn,nfft,nspden_updn,1,rhor(:,1:nspden_updn),xc_denpos)
374 
375 !There seems to be a bug in the intel11 compiler
376 !rhor = reshape(wvl_den%denspot%rhov, shape(rhor))
377  do ispden=1,nspden
378    do ifft=1,nfft
379      ii=ifft+nfft*(ispden-1)
380 !    rhor(ifft,ispden)=wvl_den%denspot%rhov(ii)
381      wvl_den%denspot%rhov(ii)=rhor(ifft,ispden)
382    end do
383  end do
384  wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY
385  write(message, '(a,a,a,a)' ) ch10, ' wvl_initro : but why are you copying me :..o('
386  call wrtout(std_out,message,'COLL')
387 
388 
389 #else
390  BIGDFT_NOTENABLED_ERROR()
391  if (.false.) write(std_out,*) me,natom,ntypat,nfft,nspden,n1,n2,n1i,n2i,n3,h(1),&
392 & pawrad(1)%mesh_size,pawtab(1)%mesh_size,spinat(1,1),zion(1),rhor(1,1),xc_denpos,&
393 & geocode,wvl_den%symObj,atindx1(1),nattyp(1),psppar(1,1,1),rprimd(1,1),xred(1,1)
394 #endif
395 
396  DBG_EXIT("COLL")
397 
398 end subroutine wvl_initro

ABINIT/wvl_mkrho [ Functions ]

[ Top ] [ Functions ]

NAME

 wvl_mkrho

FUNCTION

 This method is just a wrapper around the BigDFT routine to compute the
 density from the wavefunctions.

INPUTS

  dtset <type(dataset_type)>=input variables.
  mpi_enreg=informations about MPI parallelization
  occ(dtset%mband)=occupation numbers.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  wvl_wfs <type(wvl_projector_type)>=wavefunctions informations for wavelets.

OUTPUT

  rhor(dtset%nfft)=electron density in r space

SIDE EFFECTS

  proj <type(wvl_projector_type)>=projectors informations for wavelets.
   | proj(OUT)=computed projectors.

PARENTS

      afterscfloop,gstate,mkrho,mover,vtorho

CHILDREN

      communicate_density,sumrho,wrtout,wvl_rho_abi2big

SOURCE

431 subroutine wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den)
432 
433 #if defined HAVE_BIGDFT
434   use BigDFT_API, only : sumrho, symmetry_data, ELECTRONIC_DENSITY, communicate_density
435 #endif
436 
437 !This section has been created automatically by the script Abilint (TD).
438 !Do not modify the following lines by hand.
439 #undef ABI_FUNC
440 #define ABI_FUNC 'wvl_mkrho'
441 !End of the abilint section
442 
443  implicit none
444 
445 !Arguments -------------------------------
446 !scalars
447  type(MPI_type),intent(in) :: mpi_enreg
448  type(dataset_type),intent(in) :: dtset
449  type(wvl_wf_type),intent(inout) :: wvl_wfs
450  type(wvl_denspot_type), intent(inout) :: wvl_den
451 !arrays
452  real(dp),intent(inout) :: rhor(dtset%nfft,dtset%nspden)
453  integer, target, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,  &
454 &               (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
455  real(dp), target, intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym),  &
456 &                                 (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
457 
458 !Local variables-------------------------------
459 #if defined HAVE_BIGDFT
460 !scalars
461  character(len=500) :: message
462  integer :: comm,me,nproc
463  type(symmetry_data) :: sym
464  !for debugging:
465  !integer::ifile,ierr
466 #endif
467 
468 ! *************************************************************************
469 
470  DBG_ENTER("COLL")
471 
472 #if defined HAVE_BIGDFT
473  comm=mpi_enreg%comm_wvl
474  me=xmpi_comm_rank(comm)
475  nproc=xmpi_comm_size(comm)
476 
477  sym%symObj = wvl_den%symObj
478  sym%irrzon => irrzon
479  sym%phnons => phnons
480 
481  call sumrho(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs,wvl_wfs%ks%Lzd,&
482 & wvl_wfs%GPU,sym,wvl_den%denspot%rhod,wvl_den%denspot%xc,&
483 & wvl_wfs%ks%psi,wvl_den%denspot%rho_psi)
484 
485  call communicate_density(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs%nspin,&
486 & wvl_den%denspot%rhod,wvl_den%denspot%rho_psi,wvl_den%denspot%rhov,.false.)
487 
488  wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY
489  write(message, '(a,a,a,a)' ) ch10, ' wvl_mkrho : but why are you copying me :..o('
490  call wrtout(std_out,message,'COLL')
491 
492  call wvl_rho_abi2big(2,rhor,wvl_den)
493 
494 #else
495  BIGDFT_NOTENABLED_ERROR()
496  if (.false.) write(std_out,*) mpi_enreg%me,dtset%nstep,wvl_wfs%ks,wvl_den%symObj,&
497 & rhor(1,1),irrzon(1,1,1),phnons(1,1,1)
498 #endif
499 
500  DBG_EXIT("COLL")
501 
502 end subroutine wvl_mkrho

ABINIT/wvl_prcref [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_prcref

FUNCTION

  FIXME: add description.

INPUTS

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

SIDE EFFECTS

NOTES

PARENTS

      newrho,newvtr

CHILDREN

SOURCE

529 subroutine wvl_prcref(dielar,iprcel,my_natom,nfftprc,npawmix,nspden,pawrhoij,&
530 & rhoijrespc,usepaw,vresid,vrespc)
531 
532 
533 !This section has been created automatically by the script Abilint (TD).
534 !Do not modify the following lines by hand.
535 #undef ABI_FUNC
536 #define ABI_FUNC 'wvl_prcref'
537 !End of the abilint section
538 
539  implicit none
540 
541 !Arguments ------------------------------------
542  integer , intent(in)  :: iprcel,nfftprc,my_natom,npawmix,nspden,usepaw
543  real(dp), intent(in)  :: dielar(7)
544  real(dp), intent(in)  :: vresid(nfftprc,nspden)
545  real(dp),intent(out) :: rhoijrespc(npawmix)
546  real(dp),intent(out) :: vrespc(nfftprc,nspden)
547  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
548 
549 !Local variables-------------------------------
550  integer :: iatom,index,ispden,klmn,kmix
551  real(dp):: diemix,diemixmag,mixfac_eff
552  character(len=500) :: message                   ! to be uncommented, if needed
553 
554 ! *************************************************************************
555 
556  DBG_ENTER("COLL")
557 
558 !PENDING:
559 !optres==1, or 0, for density or potential mixing.
560 !for potential mixing, we have to average over spins.
561 !check prcref.F90 and moddiel.F90
562 
563 
564 #if defined HAVE_BIGDFT
565 #endif
566 
567  if(iprcel .ne. 0) then
568    write(message, '(a,i3,a,a,a,a)' )&
569 &   '  From the calling routine, iprcel=',iprcel,ch10,&
570 &   '  For wavelets, the only allowed value is 0.',ch10,&
571 &   '  Action : correct your input file.'
572    MSG_ERROR(message)
573  end if
574 
575 !call wvl_moddiel  !PENDING
576  diemix=dielar(4)
577 !dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ;
578  diemixmag=abs(dielar(7))
579  vrespc(:,1)=diemix*vresid(:,1)
580  if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden)
581 
582 !3) PAW only : precondition the rhoij quantities (augmentation
583 !occupancies) residuals. Use a simple preconditionning
584 !with the same mixing factor as the model dielectric function.
585 
586  if (usepaw==1.and.my_natom>0) then
587 !  mixfac=dielar(4);mixfacmag=abs(dielar(7))
588    if (pawrhoij(1)%cplex==1) then
589      index=0
590      do iatom=1,my_natom
591        do ispden=1,pawrhoij(iatom)%nspden
592          mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag
593          do kmix=1,pawrhoij(iatom)%lmnmix_sz
594            index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
595            rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
596          end do
597        end do
598      end do
599    else
600      index=-1
601      do iatom=1,my_natom
602        do ispden=1,pawrhoij(iatom)%nspden
603          mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag
604          do kmix=1,pawrhoij(iatom)%lmnmix_sz
605            index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
606            rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
607          end do
608        end do
609      end do
610    end if
611  end if
612 
613 
614 
615  DBG_EXIT("COLL")
616 
617 end subroutine wvl_prcref