TABLE OF CONTENTS


ABINIT/m_paw_mkrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_mkrho

FUNCTION

  This module contains routines used to compute PAW density on the real space fine grid.

COPYRIGHT

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

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_mkrho
24 
25  use defs_basis
26  use defs_abitypes
27  use m_abicore
28  use m_errors
29  use m_xmpi
30 
31  use m_time, only : timab
32  use m_pawang,           only : pawang_type
33  use m_pawrad,           only : pawrad_type,pawrad_deducer0
34  use m_pawtab,           only : pawtab_type,pawtab_get_lsize
35  use m_paw_sphharm,      only : initylmr
36  use m_pawfgrtab,        only : pawfgrtab_type,pawfgrtab_init,pawfgrtab_free
37  use m_pawrhoij,         only : pawrhoij_type,pawrhoij_copy,pawrhoij_free_unpacked, &
38 &                               pawrhoij_nullify,pawrhoij_free,symrhoij
39  use m_pawfgr,           only : pawfgr_type
40  use m_paw_nhat,         only : pawmknhat,nhatgrid
41  use m_paral_atom,       only : get_my_atmtab,free_my_atmtab
42  use m_fourier_interpol, only : transgrid
43 
44  use m_sort,             only : sort_dp
45  use m_splines,          only : spline,splint
46  use m_io_tools,         only : open_file
47  use m_geometry,         only : xred2xcart
48  use m_pptools,          only : printxsf
49  use m_fft,              only : fourdp
50 
51  implicit none
52 
53  private
54 
55 !public procedures.
56  public :: pawmkrho ! Build PAW electron density on fine grid, including compensation charge density
57  public :: denfgr   ! Build complete PAW electron density on fine grid, including on-site contributions
58 
59 CONTAINS  !========================================================================================

m_paw_mkrho/denfgr [ Functions ]

[ Top ] [ m_paw_mkrho ] [ Functions ]

NAME

 denfgr

FUNCTION

  Construct complete electron density on fine grid, by removing nhat
  and adding PAW corrections

INPUTS

   atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
   gmet(3,3)=reciprocal space metric tensor in bohr**-2.
   mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
   comm_atom=--optional-- MPI communicator over atoms
   my_natom=number of atoms treated by current processor
   natom= number of atoms in cell
   nattyp(ntypat)= # atoms of each type.
   ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
   nhat(pawfgr%nfft,nspden)= compensation charge density used in PAW
   nspinor=Number of spinor components
   nsppol=Number of independent spin components.
   nspden= number of spin densities
   ntypat= number of types of atoms in the cell
   pawfgr <type(pawfgr_type)>= data about the fine grid
   pawrad(ntypat) <type(pawrad_type)>= radial mesh data for each type of atom
   pawrhoij(natom) <type(pawrhoij_type)>= rho_ij data for each atom
   pawtab(ntypat) <type(pawtab_type)>= PAW functions around each type of atom
   rhor(pawfgr%nfft,nspden)= input density ($\tilde{n}+\hat{n}$ in PAW case)
   rprimd(3,3)=dimensional primitive translations for real space (bohr)
   typat(natom)= list of atom types
   ucvol=unit cell volume (bohr**3)
   xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

 rhor_paw(pawfgr%nfft,nspden)= full electron density on the fine grid

NOTES

   In PAW calculations, the valence density present in rhor includes the
   compensation charge density $\hat{n}$, and also doesn't include the on-site
   PAW contributions. For post-processing and proper visualization it is necessary
   to use the full electronic density, which is what this subroutine constructs.
   Specifically, it removes $\hat{n}$ from rhor, and also computes the on-site PAW
   terms. This is nothing other than the proper PAW treatment of the density
   operator $|\mathbf{r}\rangle\langle\mathbf{r}|$, and yields the formula
   $$\tilde{n}+\sum_{ij}\rho_ij\left[\varphi_i(\mathbf{r})\varphi_j(\mathbf{r})-
   \tilde{\varphi}_i(\mathbf{r})\tilde{\varphi}_j(\mathbf{r})\right]$$
   Notice that this formula is expressed on the fine grid, and requires
   interpolating the PAW radial functions onto this grid, as well as calling
   initylmr in order to get the angular functions on the grid points.

PARENTS

      bethe_salpeter,outscfcv,sigma

CHILDREN

      free_my_atmtab,get_my_atmtab,initylmr,nhatgrid,pawfgrtab_free
      pawfgrtab_init,pawrad_deducer0,pawtab_get_lsize,printxsf,sort_dp,spline
      splint,wrtout,xmpi_barrier,xmpi_sum,xred2xcart

SOURCE

349  subroutine denfgr(atindx1,gmet,spaceComm_in,my_natom,natom,nattyp,ngfft,nhat,nspinor,nsppol,nspden,ntypat, &
350 & pawfgr,pawrad,pawrhoij,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,rhor_nt_one,rprimd,typat,ucvol,xred,&
351 & abs_n_tilde_nt_diff,znucl,mpi_atmtab,comm_atom) ! Optional arguments
352 
353 
354 !This section has been created automatically by the script Abilint (TD).
355 !Do not modify the following lines by hand.
356 #undef ABI_FUNC
357 #define ABI_FUNC 'denfgr'
358 !End of the abilint section
359 
360  implicit none
361 
362 !Arguments ------------------------------------
363 !scalars
364  integer,intent(in) :: my_natom,natom,nspden,ntypat,prtvol,nsppol,nspinor
365  integer,optional,intent(in) :: comm_atom
366  real(dp),intent(in) :: ucvol
367  type(pawfgr_type),intent(in) :: pawfgr
368 !arrays
369  integer,intent(in) :: spaceComm_in
370  integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),typat(natom)
371  integer,optional,target,intent(in) :: mpi_atmtab(:)
372  real(dp),intent(in) :: gmet(3,3),nhat(pawfgr%nfft,nspden)
373  real(dp),intent(in) :: rhor(pawfgr%nfft,nspden),rprimd(3,3)
374  real(dp),intent(inout) :: xred(3,natom)
375  real(dp),intent(out) :: rhor_paw(pawfgr%nfft,nspden)
376  real(dp),intent(out) :: rhor_n_one(pawfgr%nfft,nspden)
377  real(dp),intent(out) :: rhor_nt_one(pawfgr%nfft,nspden)
378  real(dp),optional,intent(out) :: abs_n_tilde_nt_diff(nspden)
379  real(dp),optional,intent(in) :: znucl(ntypat)
380  type(pawrad_type),intent(in) :: pawrad(ntypat)
381  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
382  type(pawtab_type),target,intent(in) :: pawtab(ntypat)
383 
384 !Local variables-------------------------------
385 !scalars
386  integer,parameter :: master=0
387  integer :: delta,iatom,ierr,ifgd,ifftsph,inl,inrm,ipsang,irhoij
388  integer :: ispden,itypat,il,im,ilm,iln,ilmn
389  integer :: jl,jlm,jln,jm,j0lmn,jlmn
390  integer :: klmn,my_comm_atom,my_start_indx,my_end_indx
391  integer :: nfgd,nnl,normchoice,nprocs,optcut,optgr0,optgr1,optgr2
392  integer :: optrad,option,my_rank,remainder,tmp_unt
393  real(dp) :: phj,phi,rR,tphj,tphi,ybcbeg,ybcend
394  logical :: my_atmtab_allocated,paral_atom
395  character(len=500) :: message
396 !arrays
397  integer,allocatable :: l_size_atm(:),nrm_ifftsph(:)
398  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
399  integer,pointer :: my_atmtab(:)
400  real(dp) :: ylmgr(3,3,0)
401  real(dp) :: yvals(4),xcart(3,natom)
402  real(dp),allocatable :: diag(:),nrm(:),phigrd(:,:),tphigrd(:,:),ylm(:,:),ypp(:)
403  real(dp),allocatable :: phi_at_zero(:),tphi_at_zero(:)
404  real(dp),allocatable :: rhor_tmp(:,:),tot_rhor(:)
405  character(len=fnlen) :: xsf_fname
406  type(pawfgrtab_type) :: local_pawfgrtab(my_natom)
407 
408 ! ************************************************************************
409 
410  DBG_ENTER("COLL")
411 
412 !Set up parallelism over atoms (compatible only with band-FFT parallelism)
413  paral_atom=(present(comm_atom).and.(my_natom/=natom))
414  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
415  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
416  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
417 & my_natom_ref=my_natom)
418 
419 !MG: FIXME It won't work if atom-parallelism is used
420 !but even the loop over atoms below should be rewritten in this case.
421 
422 !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres
423 !the usual pawfgrtab uses r_shape which may not be the same as r_paw
424  if (my_natom>0) then
425    if (paral_atom) then
426      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat,mpi_atmtab=my_atmtab)
427      call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%cplex,l_size_atm,nspden,typat,&
428 &     mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
429    else
430      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,typat)
431      call pawfgrtab_init(local_pawfgrtab,pawrhoij(1)%cplex,l_size_atm,nspden,typat)
432    end if
433    ABI_DEALLOCATE(l_size_atm)
434  end if
435 
436 !Note: call to nhatgrid: comm_fft not used because FFT parallelism
437 !is done manually below
438  optcut = 1 ! use rpaw to construct local_pawfgrtab
439  optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally
440  optrad = 1 ! do store r-R
441  if (paral_atom) then
442    call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
443 &   optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred,&
444 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
445  else
446    call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
447 &   optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,rprimd,typat,ucvol,xred)
448  end if
449 !now local_pawfgrtab is ready to use
450 
451 !Initialise output arrays.
452  rhor_paw=zero; rhor_n_one=zero; rhor_nt_one=zero
453 
454 !Initialise and check parallell execution
455  my_rank = xmpi_comm_rank(spaceComm_in)
456  nprocs = xmpi_comm_size(spaceComm_in)
457 
458 
459 !loop over atoms in cell
460  do iatom = 1, my_natom
461    itypat = pawrhoij(iatom)%itypat
462    indlmn => pawtab(itypat)%indlmn
463    nfgd = local_pawfgrtab(iatom)%nfgd ! number of points in the fine grid for this PAW sphere
464    nnl = pawtab(itypat)%basis_size ! number of nl elements in PAW basis
465 
466 !  Division of fine grid points among processors
467    if (nprocs==1) then ! Make sure everything runs with one proc
468      write(message,'(a)') '  In denfgr - number of processors:     1'
469      call wrtout(std_out,message,'COLL')
470      write(message,'(a)') '  Calculation of PAW density done in serial'
471      call wrtout(std_out,message,'COLL')
472      write(message,'(a,I6)') '  Number of fine grid points:',nfgd
473      call wrtout(std_out,message,'COLL')
474      my_start_indx = 1
475      my_end_indx = nfgd
476    else ! Divide up the fine grid points among the processors
477      write(message,'(a,I4)') '  In denfgr - number of processors: ',nprocs
478      call wrtout(std_out,message,'COLL')
479      write(message,'(a)') '  Calculation of PAW density done in parallel'
480      call wrtout(std_out,message,'COLL')
481      write(message,'(a,I6)') '  Number of fine grid points:',nfgd
482      call wrtout(std_out,message,'COLL')
483 !    Divide the fine grid points among the processors
484      delta = int(floor(real(nfgd)/real(nprocs)))
485      remainder = nfgd-nprocs*delta
486      my_start_indx = 1+my_rank*delta
487      my_end_indx = (my_rank+1)*delta
488 !    Divide the remainder points among the processors
489 !    by shuffling indices
490      if ((my_rank+1)>remainder) then
491        my_start_indx = my_start_indx + remainder
492        my_end_indx = my_end_indx + remainder
493      else
494        my_start_indx = my_start_indx + my_rank
495        my_end_indx = my_end_indx + my_rank + 1
496      end if
497      if (prtvol>9) then
498        write(message,'(a,I6)') '  My index Starts at: ',my_start_indx
499        call wrtout(std_out,message,'PERS')
500        write(message,'(a,I6)') '             Ends at: ',my_end_indx
501        call wrtout(std_out,message,'PERS')
502        write(message,'(a,I6)') '               # pts: ',my_end_indx+1-my_start_indx
503        call wrtout(std_out,message,'PERS')
504      end if
505    end if
506 
507    write(message,'(a,I3,a,I3)') '  Entered loop for atom: ',iatom,' of:',natom
508    call wrtout(std_out,message,'PERS')
509 
510 !  obtain |r-R| values on fine grid
511    ABI_ALLOCATE(nrm,(nfgd))
512    do ifgd=1, nfgd
513      nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom)%rfgd(:,ifgd),local_pawfgrtab(iatom)%rfgd(:,ifgd)))
514    end do ! these are the |r-R| values
515 
516 !  compute Ylm for each r-R vector.
517 !  ----
518    ipsang = 1 + (pawtab(itypat)%l_size - 1)/2 ! recall l_size=2*l_max+1
519    ABI_ALLOCATE(ylm,(ipsang*ipsang,nfgd))
520    option = 1 ! compute Ylm(r-R) for vectors
521    normchoice = 1 ! use computed norms of input vectors
522    call initylmr(ipsang,normchoice,nfgd,nrm,option,local_pawfgrtab(iatom)%rfgd,ylm,ylmgr)
523 
524 !  in order to do spline fits, the |r-R| data must be sorted
525 !  ----
526    ABI_ALLOCATE(nrm_ifftsph,(nfgd))
527    nrm_ifftsph(:) = local_pawfgrtab(iatom)%ifftsph(:) ! copy of indices of points, to be rearranged by sort_dp
528    call sort_dp(nfgd,nrm,nrm_ifftsph,tol8) ! sort the nrm points, keeping track of which goes where
529 
530 !  now make spline fits of phi and tphi  onto the fine grid around the atom
531 !  ----
532    ABI_ALLOCATE(phigrd,(nfgd,nnl))
533    ABI_ALLOCATE(tphigrd,(nfgd,nnl))
534    ABI_ALLOCATE(phi_at_zero,(nnl))
535    ABI_ALLOCATE(tphi_at_zero,(nnl))
536    ABI_ALLOCATE(ypp,(pawtab(itypat)%mesh_size))
537    ABI_ALLOCATE(diag,(pawtab(itypat)%mesh_size))
538 
539    do inl = 1, nnl
540 
541 !    spline phi onto points
542      ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero;
543      call spline(pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp)
544      call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%phi(:,inl),ypp,nfgd,nrm,phigrd(:,inl))
545 
546 !    next splint tphi onto points
547      ypp(:) = zero; diag(:) = zero; ybcbeg = zero; ybcend = zero;
548      call spline(pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),pawtab(itypat)%mesh_size,ybcbeg,ybcend,ypp)
549      call splint(pawtab(itypat)%mesh_size,pawrad(itypat)%rad,pawtab(itypat)%tphi(:,inl),ypp,nfgd,nrm,tphigrd(:,inl))
550 
551 !    Find out the value of the basis function at zero using extrapolation
552      yvals = zero
553 !    Extrapolate only if this is an s-state (l=0)
554      if (indlmn(1,inl)==0) then
555        yvals(2:4) = pawtab(itypat)%phi(2:4,inl)/pawrad(itypat)%rad(2:4)
556        call pawrad_deducer0(yvals,4,pawrad(itypat))
557        write(std_out,*) 'phi_at_zero: ',yvals(1),' from:',yvals(2:4)
558      end if
559      phi_at_zero(inl) = yvals(1)
560 
561      yvals = zero
562 !    Extrapolate only if this is an s-state (l=0)
563      if (indlmn(1,inl)==0) then
564        yvals(2:4) = pawtab(itypat)%tphi(2:4,inl)/pawrad(itypat)%rad(2:4)
565        call pawrad_deducer0(yvals,4,pawrad(itypat))
566        write(std_out,*) 'tphi_at_zero: ',yvals(1),' from:',yvals(2:4)
567      end if
568      tphi_at_zero(inl) = yvals(1)
569 
570    end do ! end loop over nnl basis functions
571    ABI_DEALLOCATE(ypp)
572    ABI_DEALLOCATE(diag)
573 
574 !  loop over basis elements for this atom
575 !  because we have to store things like <phi|r'><r'|phi>-<tphi|r'><r'|tphi> at each point of the
576 !  fine grid, there is no integration, and hence no simplifications of the Y_lm's. That's why
577 !  we have to loop through the basis elements in exhaustive detail, rather than just a loop over
578 !  lmn2_size or something comparable.
579 !  ----
580    if (prtvol>9) then
581      write(message,'(a,I3)') '  Entering j-loop over basis elements for atom:',iatom
582      call wrtout(std_out,message,'PERS')
583    end if
584 
585    do jlmn=1,pawtab(itypat)%lmn_size
586 
587      if (prtvol>9) then
588        write(message,'(2(a,I3))') '  Element:',jlmn,' of:',pawtab(itypat)%lmn_size
589        call wrtout(std_out,message,'PERS')
590      end if
591 
592      jl=indlmn(1,jlmn)
593      jm=indlmn(2,jlmn)
594      jlm=indlmn(4,jlmn)
595      jln=indlmn(5,jlmn)
596      j0lmn=jlmn*(jlmn-1)/2
597 
598      if (prtvol>9) then
599        write(message,'(a,I3)') '  Entering i-loop for j:',jlmn
600        call wrtout(std_out,message,'PERS')
601      end if
602 
603      do ilmn=1,jlmn
604 
605        if (prtvol>9) then
606          write(message,'(2(a,I3))') '    Element:',ilmn,' of:',jlmn
607          call wrtout(std_out,message,'PERS')
608        end if
609 
610        il=indlmn(1,ilmn)
611        im=indlmn(2,ilmn)
612        iln=indlmn(5,ilmn)
613        ilm=indlmn(4,ilmn)
614        klmn=j0lmn+ilmn
615 
616        if (prtvol>9) then
617          write(message,'(a)') '    Entering loop over nonzero elems of rhoij'
618          call wrtout(std_out,message,'PERS')
619        end if
620 
621 !      Loop over non-zero elements of rhoij
622        do irhoij=1,pawrhoij(iatom)%nrhoijsel
623          if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn
624 
625            do ifgd=my_start_indx, my_end_indx ! loop over fine grid points in current PAW sphere
626              ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
627 
628 !            have to retrieve the spline point to use since these were sorted
629              do inrm=1, nfgd
630                if(nrm_ifftsph(inrm) == ifftsph) exit ! have found nrm point corresponding to nfgd point
631              end do ! now inrm is the index of the sorted nrm vector to use
632 
633 !            avoid division by zero
634              if(nrm(inrm) > zero) then
635                rR = nrm(inrm) ! value of |r-R| in the following
636 !              recall that <r|phi>=u(r)*Slm(r^)/r
637                phj  = phigrd(inrm,jln)*ylm(jlm,ifgd)/rR
638                phi  = phigrd(inrm,iln)*ylm(ilm,ifgd)/rR
639                tphj = tphigrd(inrm,jln)*ylm(jlm,ifgd)/rR
640                tphi = tphigrd(inrm,iln)*ylm(ilm,ifgd)/rR
641              else
642 !              use precalculated <r|phi>=u(r)*Slm(r^)/r at r=0
643                phj  = phi_at_zero(jln)*ylm(jlm,ifgd)
644                phi  = phi_at_zero(iln)*ylm(ilm,ifgd)
645                tphj = tphi_at_zero(jln)*ylm(jlm,ifgd)
646                tphi = tphi_at_zero(iln)*ylm(ilm,ifgd)
647              end if ! check if |r-R| = 0
648 
649              do ispden=1,nspden
650                if (pawrhoij(iatom)%cplex == 1) then
651                  rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + &
652 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*(phj*phi - tphj*tphi)
653 
654                  rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + &
655 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*phj*phi
656 
657                  rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + &
658 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*tphj*tphi
659                else
660                  rhor_paw(ifftsph,ispden) = rhor_paw(ifftsph,ispden) + &
661 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*(phj*phi - tphj*tphi)
662 
663                  rhor_n_one(ifftsph,ispden) = rhor_n_one(ifftsph,ispden) + &
664 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*phj*phi
665 
666                  rhor_nt_one(ifftsph,ispden) = rhor_nt_one(ifftsph,ispden) + &
667 &                 pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*tphj*tphi
668                end if ! end check on cplex rhoij
669 
670              end do ! end loop over nsdpen
671            end do ! end loop over nfgd
672          end if ! end selection on rhoij /= 0
673        end do ! end loop over non-zero rhoij
674      end do ! end loop over ilmn atomic basis states
675    end do ! end loop over jlmn atomic basis states
676 
677    ABI_DEALLOCATE(nrm)
678    ABI_DEALLOCATE(nrm_ifftsph)
679    ABI_DEALLOCATE(phigrd)
680    ABI_DEALLOCATE(tphigrd)
681    ABI_DEALLOCATE(ylm)
682    ABI_DEALLOCATE(phi_at_zero)
683    ABI_DEALLOCATE(tphi_at_zero)
684  end do     ! Loop on atoms
685 
686 !MPI sum on each node the different contributions to the PAW densities.
687  call xmpi_sum(rhor_paw,spaceComm_in,ierr)
688  call xmpi_sum(rhor_n_one,spaceComm_in,ierr)
689  call xmpi_sum(rhor_nt_one,spaceComm_in,ierr)
690  if (paral_atom) then
691    call xmpi_sum(rhor_paw,my_comm_atom,ierr)
692    call xmpi_sum(rhor_n_one,my_comm_atom,ierr)
693    call xmpi_sum(rhor_nt_one,my_comm_atom,ierr)
694  end if
695 
696  call wrtout(std_out,' *** Partial contributions to PAW rhor summed ***','PERS')
697  call xmpi_barrier(spaceComm_in)
698 
699 !Add the plane-wave contribution \tilde{n} and remove \hat{n}
700 !BE careful here since the storage mode of rhoij and rhor is different.
701  select case (nspinor)
702  case (1)
703    if (nsppol==1)  then
704      rhor_paw = rhor_paw + rhor - nhat
705    else  ! Spin-polarised case: rhor_paw contains rhor_paw(spin_up,spin_down) but we need rhor_paw(total,spin_up)
706      ABI_ALLOCATE(tot_rhor,(pawfgr%nfft))
707 !
708 !      AE rhor
709      tot_rhor(:) = SUM(rhor_paw,DIM=2)
710      rhor_paw(:,2) = rhor_paw(:,1)
711      rhor_paw(:,1) = tot_rhor
712      rhor_paw = rhor_paw + rhor - nhat
713 !
714 !      onsite AE rhor
715      tot_rhor(:) = SUM(rhor_n_one,DIM=2)
716      rhor_n_one(:,2) = rhor_n_one(:,1)
717      rhor_n_one(:,1) = tot_rhor
718 !
719 !      onsite PS rhor
720      tot_rhor(:) = SUM(rhor_nt_one,DIM=2)
721      rhor_nt_one(:,2) = rhor_nt_one(:,1)
722      rhor_nt_one(:,1) = tot_rhor
723 
724      ABI_DEALLOCATE(tot_rhor)
725    end if
726 
727  case (2)
728 !    * if nspden==4, rhor contains (n^11, n^22, Re[n^12], Im[n^12].
729 !    Storage mode for rhoij is different, See pawaccrhoij.
730    MSG_ERROR("nspinor 2 not coded")
731  case default
732    write(message,'(a,i0)')" Wrong value for nspinor=",nspinor
733    MSG_ERROR(message)
734  end select
735 
736 !if (prtvol>9) then ! Check normalisation
737 !write(message,'(a,F8.4)') '  PAWDEN - NORM OF DENSITY: ',SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
738 !call wrtout(std_out,message,'COLL')
739 !end if
740 
741  if (present(abs_n_tilde_nt_diff).AND.present(znucl)) then
742    ABI_ALLOCATE(rhor_tmp,(pawfgr%nfft,nspden))
743    do ispden=1,nspden
744      rhor_tmp(:,ispden) = zero
745      do iatom=1,my_natom
746        do ifgd=1,local_pawfgrtab(iatom)%nfgd ! loop over fine grid points in current PAW sphere
747          ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
748          rhor_tmp(ifftsph,ispden) = rhor(ifftsph,ispden) - nhat(ifftsph,ispden) &
749 &         - rhor_nt_one(ifftsph,ispden)
750        end do !ifgd
751      end do ! iatom
752    end do ! ispden
753    if (paral_atom) then
754      call xmpi_sum(rhor_tmp,my_comm_atom,ierr)
755    end if
756 
757    if (my_rank==master) then
758      do ispden=1,nspden
759 !      Write to xsf file
760        call xred2xcart(natom,rprimd,xcart,xred)
761        write(xsf_fname,'(a,I0,a)') 'N_tilde_onsite_diff_sp',ispden,'.xsf'
762        if (open_file(xsf_fname,message, unit=tmp_unt,status='unknown',form='formatted') /= 0) then
763          MSG_ERROR(message)
764        end if
765        call printxsf(ngfft(1),ngfft(2),ngfft(3),rhor_tmp(:,ispden),rprimd,&
766 &       (/zero,zero,zero/),natom,ntypat,typat,xcart,znucl,tmp_unt,0)
767        close(tmp_unt)
768        abs_n_tilde_nt_diff(ispden) = SUM(ABS(rhor_tmp(:,ispden)))/pawfgr%nfft
769        write(message,'(4(a),F16.9,2(a,I0),a)') ch10,'  Wrote xsf file with \tilde{n}-\tilde{n}^1.',ch10,&
770 &       '  Value of norm |\tilde{n}-\tilde{n}^1|:',&
771 &       abs_n_tilde_nt_diff(ispden),' spin: ',ispden,' of ',nspden,ch10
772        call wrtout(std_out,message,'COLL')
773      end do
774    end if
775    ABI_DEALLOCATE(rhor_tmp)
776 
777  else if ((present(abs_n_tilde_nt_diff).AND.(.NOT.present(znucl))) &
778 &   .OR.(.NOT.present(abs_n_tilde_nt_diff).AND.(present(znucl)))) then
779    write(message,'(a)') ' Both abs_n_tilde_nt_diff *and* znucl must be passed',ch10,&
780 &   'to denfgr for |\tilde{n}-\tilde{n}^1| norm evaluation.'
781    MSG_ERROR(message)
782  end if
783 
784  call pawfgrtab_free(local_pawfgrtab)
785 
786 !Destroy atom table used for parallelism
787  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
788 
789  call xmpi_barrier(spaceComm_in)
790 
791  DBG_EXIT("COLL")
792 
793  end subroutine denfgr

m_paw_mkrho/pawmkrho [ Functions ]

[ Top ] [ m_paw_mkrho ] [ Functions ]

NAME

 pawmkrho

FUNCTION

 PAW only:
 Build total pseudo (compensated) density (\tild_rho + \hat_rho)
 Build compensation charge density (\hat_rho)
 Build occupation matrix (packed storage)

INPUTS

  compute_rhor_rhog: if 1: set the computation of rhor and rhog in addition to the compensating charge.
                     if 0: compute only the compensating charge
  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
         1 for GS calculations
  gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1).
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  ipert=index of perturbation if pawrhoij is a pertubed rhoij
        no meaning for ground-state calculations (should be 0)
  idir=direction of atomic displacement (in case of atomic displ. perturb.)
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell
  nspden=number of spin-density components
  nsym=number of symmetry elements in space group
  ntypat=number of types of atoms in unit cell.
  paral_kgb=option for (kpt,g vectors,bands) parallelism
  pawang <type(pawang_type)>=angular mesh discretization and related data
  pawang_sym <type(pawang_type)>=angular data used for symmetrization
                                 optional parameter only needed for RF calculations
  pawfgr <type(paw_fgr_type)>=fine rectangular grid parameters
  pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawprtvol=control print volume and debugging output for PAW
  pawrhoij0(natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0)
                                          optional parameter only needed for RF calculations
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  qphon(3)=wavevector of the phonon (RF only)
  rhopsg(2,pawfgr%nfftc)= pseudo density given on the coarse grid in reciprocal space
  rhopsr(pawfgr%nfftc,nspden)= pseudo density given on the coarse grid in real space
  rprimd(3,3)=real space primitive translations.
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symrec(3,3,nsym)=symmetries of group in terms of operations on
                   reciprocal space primitive translations
  typat(natom)=type for each atom
  ucvol=volume of the unit cell
  xred(3,natom)= reduced atomic coordinates

OUTPUT

  compch_fft=compensation charge inside spheres integrated over fine fft grid
  pawnhat(pawfgr%nfft,nspden)=compensation charge density on fine rectangular grid (optional argument)
  rhog(2,pawfgr%nfft)= compensated pseudo density given on the fine grid in reciprocal space
                       This output is optional
  rhor(pawfgr%nfft,nspden)= compensated pseudo density given on the fine grid in real space

SIDE EFFECTS

  pawrhoij(my_natom)= PAW occupancies
                   At input : values at previous step  in packed storage (pawrhoij()%rhoijp)
                   At output: values (symmetrized)     in packed storage (pawrhoij()%rhoijp)
  pawrhoij_unsym(:)= unsymmetrized PAW occupancies
                   At input : values (unsymmetrized) in unpacked storage (pawrhoij()%rhoij_)
                   At output: values in unpacked storage (pawrhoij()%rhoij_) are destroyed

NOTES

  pawrhoij and pawrhoij_unsym can be identical (refer to the same pawrhoij datastructure).
  They should be different only if pawrhoij is distributed over atomic sites
  (in that case pawrhoij_unsym should not be distributed over atomic sites).

PARENTS

      afterscfloop,dfpt_nstpaw,dfpt_rhofermi,dfpt_vtorho,scfcv,vtorho

CHILDREN

      fourdp,pawmknhat,pawrhoij_copy,pawrhoij_free,pawrhoij_free_unpacked
      pawrhoij_nullify,symrhoij,timab,transgrid

SOURCE

140 subroutine pawmkrho(compute_rhor_rhog,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
141 &          my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawang,pawfgr,pawfgrtab,pawprtvol,&
142 &          pawrhoij,pawrhoij_unsym,&
143 &          pawtab,qphon,rhopsg,rhopsr,rhor,rprimd,symafm,symrec,typat,ucvol,usewvl,xred,&
144 &          pawang_sym,pawnhat,pawrhoij0,rhog) ! optional arguments
145 
146 
147 !This section has been created automatically by the script Abilint (TD).
148 !Do not modify the following lines by hand.
149 #undef ABI_FUNC
150 #define ABI_FUNC 'pawmkrho'
151 !End of the abilint section
152 
153  implicit none
154 
155 !Arguments ------------------------------------
156 !scalars
157  integer,intent(in) :: compute_rhor_rhog,cplex,idir,ipert,my_natom,natom,nspden,nsym,ntypat,paral_kgb,pawprtvol
158  integer,intent(in) :: usewvl
159  real(dp),intent(in) :: ucvol
160  real(dp),intent(out) :: compch_fft
161  type(MPI_type),intent(in) :: mpi_enreg
162  type(pawang_type),intent(in) :: pawang
163  type(pawang_type),intent(in),optional :: pawang_sym
164  type(pawfgr_type),intent(in) :: pawfgr
165 !arrays
166  integer,intent(in) :: indsym(4,nsym,natom)
167  integer,intent(in) :: symafm(nsym),symrec(3,3,nsym),typat(natom)
168  real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom)
169  real(dp),intent(inout),target,optional :: pawnhat(cplex*pawfgr%nfft,nspden) !vz_i
170  real(dp),intent(inout) :: rhor(cplex*pawfgr%nfft,nspden*compute_rhor_rhog)
171  real(dp),intent(out),optional :: rhog(2,pawfgr%nfft*compute_rhor_rhog)
172  real(dp),intent(inout) :: rhopsg(2,pawfgr%nfftc*compute_rhor_rhog),rhopsr(cplex*pawfgr%nfftc,nspden*compute_rhor_rhog)
173  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
174  type(pawrhoij_type),intent(inout),target :: pawrhoij(:)
175  type(pawrhoij_type),intent(inout) :: pawrhoij_unsym(:)
176  type(pawrhoij_type),intent(in),target,optional :: pawrhoij0(my_natom)
177  type(pawtab_type),intent(in) :: pawtab(ntypat)
178 
179 !Local variables-------------------------------
180 !scalars
181  integer :: choice,ider,izero,option
182  character(len=500) :: msg
183 !arrays
184  real(dp) :: tsec(2)
185  real(dp) :: rhodum(0,0,0)
186  real(dp),pointer :: pawnhat_ptr(:,:)
187  type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij0_ptr(:)
188 
189 ! ***********************************************************************
190 
191  DBG_ENTER("COLL")
192 
193  call timab(556,1,tsec)
194 
195 !Compatibility tests
196  if (size(pawrhoij_unsym)>0) then
197    if (pawrhoij_unsym(1)%use_rhoij_==0) then
198      msg='  rhoij_ field must be allocated in pawrhoij_unsym !'
199      MSG_BUG(msg)
200    end if
201  end if
202  if (ipert>0.and.(.not.present(pawrhoij0))) then
203    msg='  pawrhoij0 must be present when ipert>0 !'
204    MSG_BUG(msg)
205  end if
206  
207 !Symetrize PAW occupation matrix and store it in packed storage
208  call timab(557,1,tsec)
209  option=1;choice=1
210  if (present(pawang_sym)) then
211    call symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,natom,nsym,ntypat,&
212 &   option,pawang_sym,pawprtvol,pawtab,rprimd,symafm,symrec,typat,&
213 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
214 &   qphon=qphon)
215  else
216    call symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,ipert,natom,nsym,ntypat,&
217 &   option,pawang,pawprtvol,pawtab,rprimd,symafm,symrec,typat,&
218 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
219 &   qphon=qphon)
220  end if
221  call pawrhoij_free_unpacked(pawrhoij_unsym)
222  call timab(557,2,tsec)
223 
224 !In somes cases (parallelism), has to distribute the PAW occupation matrix
225  if (size(pawrhoij)==natom.and.(my_natom/=natom)) then
226    ABI_DATATYPE_ALLOCATE(pawrhoij_ptr,(my_natom))
227    call pawrhoij_nullify(pawrhoij_ptr)
228    call pawrhoij_copy(pawrhoij,pawrhoij_ptr,&
229 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom, &
230 &   keep_cplex=.false.,keep_itypat=.false.,keep_nspden=.false.)       
231  else
232    pawrhoij_ptr=>pawrhoij
233  end if
234 
235 !Compute compensation charge density
236  ider=0;izero=0
237  if (present(pawnhat)) then
238    pawnhat_ptr => pawnhat
239  else
240    ABI_ALLOCATE(pawnhat_ptr,(pawfgr%nfft,nspden))
241  end if
242  if (present(pawrhoij0)) then
243    pawrhoij0_ptr => pawrhoij0
244  else
245    pawrhoij0_ptr => pawrhoij_ptr
246  end if
247 
248  call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
249 & pawfgr%nfft,pawfgr%ngfft,ider,nspden,ntypat,pawang,pawfgrtab,&
250 & rhodum,pawnhat_ptr,pawrhoij_ptr,pawrhoij0_ptr,pawtab,qphon,rprimd,ucvol,usewvl,xred,&
251 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
252 & comm_fft=mpi_enreg%comm_fft,paral_kgb=paral_kgb,me_g0=mpi_enreg%me_g0,&
253 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 
254 
255  if (compute_rhor_rhog/=0) then
256 !  Transfer pseudo density from coarse grid to fine grid
257    if(usewvl==0) then
258      call transgrid(cplex,mpi_enreg,nspden,+1,1,0,paral_kgb,pawfgr,rhopsg,rhodum,rhopsr,rhor)
259    end if
260 
261 !  Add pseudo density and compensation charge density (on fine grid)
262    rhor(:,:)=rhor(:,:)+pawnhat_ptr(:,:)
263 
264 !  Compute compensated pseudo density in reciprocal space
265    if (present(rhog)) then
266      call fourdp(cplex,rhog,rhor(:,1),-1,mpi_enreg,pawfgr%nfft,pawfgr%ngfft,paral_kgb,0)
267    end if
268  end if
269 
270 !Free temporary memory spaces
271  if (.not.present(pawnhat)) then
272    ABI_DEALLOCATE(pawnhat_ptr)
273  end if
274  if (size(pawrhoij)==natom.and.(my_natom/=natom)) then
275    call pawrhoij_free(pawrhoij_ptr)
276    ABI_DATATYPE_DEALLOCATE(pawrhoij_ptr)
277  end if
278  nullify(pawnhat_ptr)
279  nullify(pawrhoij_ptr)
280 
281  call timab(556,2,tsec)
282 
283  DBG_EXIT("COLL")
284 
285 end subroutine pawmkrho