TABLE OF CONTENTS


ABINIT/dfpt_accrho [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_accrho

FUNCTION

 Response function calculation only:
  Accumulate contribution to first-order density due to current (k,band)
  Also accumulate zero-order potential part of the 2nd-order total energy (if needed)

INPUTS

  counter=counter for status file
  cplex=1 if 1st-order density is real, 2 if 1st-order density is complex
  cwave0(2,npw*nspinor)=GS wavefunction at k, in reciprocal space
  cwave1(2,npw1*nspinor)=1st-order wavefunction at k,q, in reciprocal space
  cwavef(2,npw1*nspinor)=1st-order wavefunction at k,q, in reciprocal space, without correction due to occupation change
  cwaveprj0(natom,nspinor*usecprj)= GS wave function at k projected with nl projectors
  cwaveprj1(natom,nspinor*usecprj)= 1st-order wave function at k,q projected with nl projectors
  filstat=name of the status file
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  iband=index of current band
  idir=direction of the current perturbation
  ipert=type of the perturbation
  isppol=1 index of current spin component
  kptopt=option for the generation of k points
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nband_k=number of bands at this k point for that spin polarization
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  npw_k=number of planewaves in basis sphere at k
  npw1_k=number of planewaves in basis sphere at k+q
  nspinor=number of spinorial components of the wavefunctions
  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
  option= 1: accumulate 1st-order density,
          2: accumulate 0-order potential part of the 2nd-order total energy
          3: accumulate both
  prtvol=control print volume and debugging output
  tim_fourwf= timing code for fourwf (5 from dfpt_vtowfk, 18 from dfpt_nstwf)
  wf_corrected=flag put to 1 if cwave1 is different from cwavef (if there is a contribution from occ. change)
  wtk_k=weight assigned to the k point.

OUTPUT

  ====== if option=2 or option=3 =====
  eloc0_k=zero-order local contribution to 2nd-order total energy for current band and k

SIDE EFFECTS

  ====== if option=1 or option=3 =====
    rhoaug1(cplex*n4,n5,n6,nvloc)= density in electrons/bohr**3,
    ==== if gs_hamkq%usepaw=1 =====
    pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
                                            (cumulative, so input as well as output)

NOTES

  In this part of the treatment of one band, one has to
  perform Fourier transforms, and to treat separately the
  two spinorial components of the wavefunction.
  Was part of dfpt_vtowfk before.

PARENTS

      dfpt_nstpaw,dfpt_vtowfk,dfpt_wfkfermi

CHILDREN

      fourwf,get_my_atmtab,getcprj,pawaccrhoij,pawcprj_alloc,pawcprj_free
      status

SOURCE

 615 subroutine dfpt_accrho(counter,cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,&
 616 &                  eloc0_k,filstat,gs_hamkq,iband,idir,ipert,isppol,kptopt,&
 617 &                  mpi_enreg,natom,nband_k,ncpgr,npw_k,npw1_k,nspinor,occ_k,&
 618 &                  option,pawrhoij1,prtvol,rhoaug1,tim_fourwf,wf_corrected,&
 619 &                  wtk_k,comm_atom,mpi_atmtab)
 620 
 621 
 622 !This section has been created automatically by the script Abilint (TD).
 623 !Do not modify the following lines by hand.
 624 #undef ABI_FUNC
 625 #define ABI_FUNC 'dfpt_accrho'
 626 !End of the abilint section
 627 
 628  implicit none
 629 
 630 !Arguments ------------------------------------
 631 !scalars
 632  integer,intent(in) :: counter,cplex,iband,idir,ipert,isppol,kptopt,natom,nband_k
 633  integer,intent(in) :: ncpgr,npw_k,npw1_k,nspinor,option,prtvol,tim_fourwf,wf_corrected
 634  integer,optional,intent(in) :: comm_atom
 635  integer,optional,target,intent(in) :: mpi_atmtab(:)
 636  real(dp),intent(in) :: wtk_k
 637  real(dp),intent(out) :: eloc0_k
 638  character(len=fnlen),intent(in) :: filstat
 639  type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
 640  type(MPI_type),intent(in) :: mpi_enreg
 641 !arrays
 642  real(dp),intent(in),target :: cwave0(2,npw_k*nspinor),cwave1(2,npw1_k*nspinor),cwavef(2,npw1_k*nspinor)
 643  real(dp),intent(in) :: occ_k(nband_k)
 644  real(dp),intent(inout) :: rhoaug1(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc)
 645  type(pawcprj_type),intent(in) :: cwaveprj0(natom,nspinor*gs_hamkq%usecprj)
 646  type(pawcprj_type),intent(in) :: cwaveprj1(natom,nspinor*gs_hamkq%usepaw)
 647  type(pawrhoij_type),intent(inout) :: pawrhoij1(:)
 648 
 649 !Local variables-------------------------------
 650 !scalars
 651  integer,parameter :: level=14
 652  integer :: choice,cplex_cprj,i1,i2,i3,iexit,ispinor,my_comm_atom,my_natom,n1,n2,n3,option_rhoij
 653  logical :: my_atmtab_allocated,paral_atom
 654  logical :: usetimerev
 655  real(dp) :: im0,im1,re0,re1,valuer,diag,offdiag,weight
 656  real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down
 657 !arrays
 658  integer,pointer :: my_atmtab(:)
 659  real(dp) :: dummy(2,1)
 660  real(dp),allocatable :: rhoaug(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:)
 661  real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:)
 662  real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:)
 663  real(dp),pointer :: cwavef_sp(:,:),cwavef_up(:,:),cwavef_down(:,:)
 664  real(dp),pointer :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:)
 665  real(dp),pointer :: vlocal(:,:,:,:)=>null()
 666  type(pawcprj_type),allocatable :: cwaveprj_tmp(:,:)
 667 !TEST
 668 !  real(dp),save :: v1,v2
 669 !  real(dp) :: r1,r2
 670 ! *********************************************************************
 671  DBG_ENTER("COLL")
 672 
 673  if (option/=1.and.option/=2.and.option/=3) return
 674 
 675 !Initializations
 676  ABI_ALLOCATE(rhoaug,(gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc))
 677  n1=gs_hamkq%ngfft(1);n2=gs_hamkq%ngfft(2);n3=gs_hamkq%ngfft(3)
 678  if (option==2.or.option==3) eloc0_k=zero
 679  if (option==2.or.option==3) vlocal => gs_hamkq%vlocal
 680 
 681 !Loop on spinorial components
 682 ! TODO : double loop on spinors for full rhoaug1 matrix if nspden =4
 683  if (gs_hamkq%nvloc/=4) then  ! see later EB FR
 684    ABI_ALLOCATE(wfraug1,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 685 
 686    do ispinor=1,nspinor
 687 
 688      if (prtvol>=10) then
 689        call status(counter,filstat,iexit,level,'density update')
 690      end if
 691 
 692 !  Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy
 693 !  --------------------------------------------------------------------------------------------
 694 
 695 !  Fourier transform of cwavef. Here, rhoaug is a dummy variable.
 696      if (wf_corrected==0.or.option==2.or.option==3) then
 697        if (ispinor==1) then
 698          cwavef_sp => cwavef(:,1:npw1_k)
 699        else
 700          cwavef_sp => cwavef(:,1+npw1_k:2*npw1_k)
 701        end if
 702        !make an inverse FFT from cwavef_sp to wfraug1
 703        call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 704 &       gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 705 &       gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 706 &       weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 707        nullify(cwavef_sp)
 708      end if
 709 
 710 !  Compute contribution of this band to zero-order potential part of the 2nd-order total energy
 711 !  NB: this is spinor diagonal
 712      if (option==2.or.option==3) then
 713        valuer=zero
 714        do i3=1,n3
 715          do i2=1,n2
 716            do i1=1,n1
 717              valuer=valuer+vlocal(i1,i2,i3,1)*(wfraug1(1,i1,i2,i3)**2+wfraug1(2,i1,i2,i3)**2)
 718            end do
 719          end do
 720        end do
 721 
 722 !    Local potential energy of this band
 723        eloc0_k=eloc0_k+two*valuer/dble(gs_hamkq%nfft)
 724      end if ! option
 725 
 726 !  Part devoted to the accumulation of the 1st-order density
 727 !  ---------------------------------------------------------
 728      if (option==1.or.option==3) then
 729 
 730 !    Compute 1st-order WF in real space
 731 !    One needs the Fourier transform of cwave1. However, only the one of
 732 !    cwavef is available. If cwavef and cwave1 differs, this Fourier
 733 !    transform must be computed. In both case the result is in wfraug1.
 734        if (wf_corrected==1) then
 735          if (ispinor==1) then
 736            cwavef_sp => cwave1(:,1:npw1_k)
 737          else
 738            cwavef_sp => cwave1(:,1+npw1_k:2*npw1_k)
 739          end if
 740          call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 741 &         gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 742 &         gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 743 &         weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 744          nullify(cwavef_sp)
 745        end if
 746 
 747 !    Compute 0-order WF in real space
 748 ! TODO: add loop over ispinor_prime here
 749        ABI_ALLOCATE(wfraug,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 750        if (ispinor==1) then
 751          cwavef_sp => cwave0(:,1:npw_k)
 752        else
 753          cwavef_sp => cwave0(:,1+npw_k:2*npw_k)
 754        end if
 755        call fourwf(1,rhoaug,cwavef_sp,dummy,wfraug,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
 756 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 757 &       gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 758 &       weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 759        nullify(cwavef_sp)
 760 
 761 !    The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
 762        weight=two*occ_k(iband)*wtk_k/gs_hamkq%ucvol
 763 !    Accumulate 1st-order density
 764        if (cplex==2) then
 765          do i3=1,n3
 766            do i2=1,n2
 767              do i1=1,n1
 768                re0=wfraug(1,i1,i2,i3)  ; im0=wfraug(2,i1,i2,i3)
 769                re1=wfraug1(1,i1,i2,i3) ; im1=wfraug1(2,i1,i2,i3)
 770 ! TODO: check which terms (ispinor ispinorp) enter a given element of rhoaug1
 771                rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1)
 772                rhoaug1(2*i1  ,i2,i3,1)=rhoaug1(2*i1  ,i2,i3,1)+weight*(re0*im1-im0*re1)
 773              end do
 774            end do
 775          end do
 776        else
 777          do i3=1,n3
 778            do i2=1,n2
 779              do i1=1,n1
 780                rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1) &
 781 &               +weight*(wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) &
 782 &               +wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3))
 783              end do
 784            end do
 785          end do
 786        end if
 787        ABI_DEALLOCATE(wfraug)
 788      end if ! option
 789 
 790    end do ! Loop on spinorial components if nspden=1 or 2
 791 
 792    ABI_DEALLOCATE(wfraug1)
 793  else ! nvloc = 4
 794 ! The same lines of code are in 72_response/dfpt_mkrho.F90
 795 ! TODO merge these lines in a single routine??!!
 796    if (prtvol>=10) then
 797      call status(counter,filstat,iexit,level,'density update')
 798    end if
 799    ABI_ALLOCATE(wfraug1_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 800    ABI_ALLOCATE(wfraug1_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 801 
 802    wfraug1_up(:,:,:,:)=zero
 803    wfraug1_down(:,:,:,:)=zero
 804 
 805 !  Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy
 806 !  --------------------------------------------------------------------------------------------
 807 
 808 !  Fourier transform of cwavef. Here, rhoaug is a dummy variable.
 809    if (wf_corrected==0.or.option==2.or.option==3) then
 810      cwavef_up => cwavef(:,1:npw1_k) ! wfs up spin-polarized
 811      call fourwf(cplex,rhoaug(:,:,:,1),cwavef_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 812 &     gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 813 &     gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 814 &     weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 815      nullify(cwavef_up)
 816 
 817      cwavef_down => cwavef(:,1+npw1_k:2*npw1_k) ! wfs down spin-polarized
 818      call fourwf(cplex,rhoaug(:,:,:,1),cwavef_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 819 &     gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 820 &     gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 821 &     weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 822      nullify(cwavef_down)
 823    end if
 824    if (option==2.or.option==3) then
 825      valuer=zero
 826      diag=zero
 827      offdiag=zero
 828    ! EB FR 2nd term in Eq. 91 PRB52,1096 [[cite:Gonze1995]] for non-collinear magnetism
 829      do i3=1,n3
 830        do i2=1,n2
 831          do i1=1,n1
 832 !TEST DFPT
 833            diag=vlocal(i1,i2,i3,1)*(wfraug1_up(1,i1,i2,i3)**2+wfraug1_up(2,i1,i2,i3)**2)&
 834 &           +vlocal(i1,i2,i3,2)*(wfraug1_down(1,i1,i2,i3)**2+wfraug1_down(2,i1,i2,i3)**2)
 835            offdiag=(two*vlocal(i1,i2,i3,3)*((wfraug1_up(1,i1,i2,i3)*wfraug1_down(1,i1,i2,i3))+&
 836 &           (wfraug1_up(2,i1,i2,i3)*wfraug1_down(2,i1,i2,i3))))+&
 837 &           (two*vlocal(i1,i2,i3,4)*((-wfraug1_down(2,i1,i2,i3)*wfraug1_up(1,i1,i2,i3))+&
 838 &           (wfraug1_down(1,i1,i2,i3)*wfraug1_up(2,i1,i2,i3))))
 839            valuer=valuer+diag+offdiag
 840          end do
 841        end do
 842      end do
 843 !    Local potential energy of this band
 844      eloc0_k=eloc0_k+two*valuer/dble(gs_hamkq%nfft)
 845    end if ! option
 846 
 847 !  Part devoted to the accumulation of the 1st-order density
 848 !  ---------------------------------------------------------
 849 
 850 ! first order
 851 !
 852    if (option==1.or.option==3) then
 853 
 854      !SPr: condition on wf_corrected not to do FFTs of the same Bloch functions
 855      if (wf_corrected==1) then
 856        cwave1_up => cwave1(:,1:npw1_k)
 857        cwave1_down => cwave1(:,1+npw1_k:2*npw1_k)
 858        wfraug1_up(:,:,:,:)=zero
 859        wfraug1_down(:,:,:,:)=zero
 860 
 861        call fourwf(cplex,rhoaug(:,:,:,1),cwave1_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 862 &       gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 863 &       gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 864 &       weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 865        nullify(cwave1_up)
 866 
 867        call fourwf(cplex,rhoaug(:,:,:,1),cwave1_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
 868 &       gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 869 &       gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 870 &       weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 871        nullify(cwave1_down)
 872      end if
 873 
 874 
 875 ! EB FR build spinorial wavefunctions
 876 ! zero order
 877      cwave0_up => cwave0(:,1:npw_k)
 878      cwave0_down => cwave0(:,1+npw_k:2*npw_k)
 879      ABI_ALLOCATE(wfraug_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 880      ABI_ALLOCATE(wfraug_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
 881      wfraug_up(:,:,:,:)=zero
 882      wfraug_down(:,:,:,:)=zero
 883 !
 884      !density components
 885      !GS wfk Fourrier Tranform
 886      ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument
 887      call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
 888 &     gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 889 &     gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 890 &     weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 891      nullify(cwave0_up)
 892      call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
 893 &     gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
 894 &     gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,mpi_enreg%paral_kgb,tim_fourwf,&
 895 &     weight,weight,use_gpu_cuda=gs_hamkq%use_gpu_cuda)
 896      nullify(cwave0_down)
 897 !    Accumulate 1st-order density (x component)
 898      re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero
 899      re1_down=zero;im1_down=zero
 900 !    The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
 901 !    SPr: the following treatment with factor=2 is ok for perturbations not breaking the
 902 !         time reversal symmetry of the Hamiltonian (due to Kramer's degeneracy) hence
 903 !         not applicable for magnetic field perturbation (for phonons with SOC, H^(0) has
 904 !         time reversal symmetry though). The formulas below are rectified in dfpt_scfcv
 905 !         in case of broken time-reversal upon reconstructing rhor1_pq and rhor1_mq.
 906      weight=two*occ_k(iband)*wtk_k/gs_hamkq%ucvol
 907      if (cplex==2) then
 908        do i3=1,n3
 909          do i2=1,n2
 910            do i1=1,n1
 911              re0_up=wfraug_up(1,i1,i2,i3)  ;     im0_up=wfraug_up(2,i1,i2,i3)
 912              re1_up=wfraug1_up(1,i1,i2,i3) ;     im1_up=wfraug1_up(2,i1,i2,i3)
 913              re0_down=wfraug_down(1,i1,i2,i3)  ; im0_down=wfraug_down(2,i1,i2,i3)
 914              re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
 915              !SPr: in case of +q/-q calculation, the factor will be corrected later from dfpt_scfcv level
 916              !     along with the reconstruction of correct rhor1_{+q} and rhor1_{-q}
 917              !     here, rhoaug1_{sigma,sigma'} = \sum_{n,k} u1_{sigma} u0*_{sigma'} independent of the sign of q
 918              rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) !n_upup
 919              rhoaug1(2*i1  ,i2,i3,1)=rhoaug1(2*i1  ,i2,i3,1)+weight*(re0_up*im1_up-im0_up*re1_up)
 920              rhoaug1(2*i1-1,i2,i3,4)=rhoaug1(2*i1-1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
 921              rhoaug1(2*i1  ,i2,i3,4)=rhoaug1(2*i1  ,i2,i3,4)+weight*(re0_down*im1_down-im0_down*re1_down)
 922 
 923              rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+im1_up*im0_down)& !Re[m1x]
 924 &            +weight*(re1_down*re0_up+im1_down*im0_up)
 925              rhoaug1(2*i1  ,i2,i3,2)=rhoaug1(2*i1  ,i2,i3,2)+weight*(-re1_up*im0_down+im1_up*re0_down)& !Im[m1x]
 926 &            +weight*(-re1_down*im0_up+im1_down*re0_up)
 927 
 928              rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(+re1_up*im0_down-im1_up*re0_down)& !Re[m1y]
 929 &            +weight*(-re1_down*im0_up+im1_down*re0_up)
 930              rhoaug1(2*i1  ,i2,i3,3)=rhoaug1(2*i1  ,i2,i3,3)+weight*(+re1_up*re0_down+im1_up*im0_down)& !Im[m1y]
 931 &            +weight*(-re1_down*re0_up-im1_down*im0_up)
 932            end do
 933          end do
 934        end do
 935      else !cplex
 936        re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero
 937        re1_down=zero;im1_down=zero
 938        do i3=1,n3
 939          do i2=1,n2
 940            do i1=1,n1
 941              re0_up=wfraug_up(1,i1,i2,i3)  ;     im0_up=wfraug_up(2,i1,i2,i3)
 942              re1_up=wfraug1_up(1,i1,i2,i3) ;     im1_up=wfraug1_up(2,i1,i2,i3)
 943              re0_down=wfraug_down(1,i1,i2,i3)  ; im0_down=wfraug_down(2,i1,i2,i3)
 944              re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
 945 
 946              rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup
 947              rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
 948              rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down &
 949 &             +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight
 950              rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down &
 951 &             +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight
 952            end do
 953          end do
 954        end do
 955      end if !cplex
 956 
 957      ABI_DEALLOCATE(wfraug_up)
 958      ABI_DEALLOCATE(wfraug_down)
 959 
 960    end if ! option
 961 
 962    ABI_DEALLOCATE(wfraug1_up)
 963    ABI_DEALLOCATE(wfraug1_down)
 964 
 965  end if ! nvloc /= 4
 966 
 967  ABI_DEALLOCATE(rhoaug)
 968 
 969 !Part devoted to the accumulation of the 1st-order occupation matrix in PAW case
 970 ! TODO: parse for more nspden 4 dependencies on spinors
 971 ! EB FR CHECK: to be modified for non-collinear?????
 972 !-------------------------------------------------------------------------------
 973 
 974  if ((option==1.or.option==3).and.gs_hamkq%usepaw==1) then
 975 
 976 !  Set up parallelism over atoms
 977    my_natom=natom; if(gs_hamkq%usepaw==1) my_natom=size(pawrhoij1)
 978    paral_atom=(present(comm_atom).and.(my_natom/=natom))
 979    my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 980    nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 981    call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
 982 
 983    cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1
 984    option_rhoij=2;usetimerev=(kptopt>0.and.kptopt<3)
 985 
 986    if (gs_hamkq%usecprj==1) then
 987 !TEST
 988 ! if (iband==1) then
 989 !  if (nspinor==1) then
 990 !   if (iband==1) write(100+iband,*) "===================================="
 991 !   if (iband==1) write(200+iband,*) "===================================="
 992 !   if (iband==1) write(300+iband,*) "===================================="
 993 !   if (iband==1) write(400+iband,*) "===================================="
 994 !   write(100+iband,*) cwaveprj0(1,1)%cp(1,1)**2+cwaveprj0(1,1)%cp(2,1)**2,&
 995 ! &                    cwaveprj0(1,1)%cp(1,2)**2+cwaveprj0(1,1)%cp(2,2)**2,&
 996 ! &                    cwaveprj0(1,1)%cp(1,3)**2+cwaveprj0(1,1)%cp(2,3)**2,&
 997 ! &                    cwaveprj0(1,1)%cp(1,4)**2+cwaveprj0(1,1)%cp(2,4)**2,&
 998 ! &                    cwaveprj0(1,1)%cp(1,5)**2+cwaveprj0(1,1)%cp(2,5)**2,&
 999 ! &                    cwaveprj0(1,1)%cp(1,6)**2+cwaveprj0(1,1)%cp(2,6)**2,&
1000 ! &                    cwaveprj0(1,1)%cp(1,7)**2+cwaveprj0(1,1)%cp(2,7)**2,&
1001 ! &                    cwaveprj0(1,1)%cp(1,8)**2+cwaveprj0(1,1)%cp(2,8)**2
1002 !   write(200+iband,*) cwaveprj1(1,1)%cp(1,1)**2+cwaveprj1(1,1)%cp(2,1)**2,&
1003 ! &                    cwaveprj1(1,1)%cp(1,2)**2+cwaveprj1(1,1)%cp(2,2)**2,&
1004 ! &                    cwaveprj1(1,1)%cp(1,3)**2+cwaveprj1(1,1)%cp(2,3)**2,&
1005 ! &                    cwaveprj1(1,1)%cp(1,4)**2+cwaveprj1(1,1)%cp(2,4)**2,&
1006 ! &                    cwaveprj1(1,1)%cp(1,5)**2+cwaveprj1(1,1)%cp(2,5)**2,&
1007 ! &                    cwaveprj1(1,1)%cp(1,6)**2+cwaveprj1(1,1)%cp(2,6)**2,&
1008 ! &                    cwaveprj1(1,1)%cp(1,7)**2+cwaveprj1(1,1)%cp(2,7)**2,&
1009 ! &                    cwaveprj1(1,1)%cp(1,8)**2+cwaveprj1(1,1)%cp(2,8)**2
1010 !   write(300+iband,*) cwave0(1,1)**2+cwave0(2,1)**2,&
1011 !                      cwave0(1,5)**2+cwave0(2,5)**2,&
1012 !                      cwave0(1,9)**2+cwave0(2,9)**2
1013 !   write(400+iband,*) cwave1(1,1)**2+cwave1(2,1)**2,&
1014 !                      cwave1(1,5)**2+cwave1(2,5)**2,&
1015 !                      cwave1(1,9)**2+cwave1(2,9)**2
1016 !  end if
1017 !  if (nspinor==2) then
1018 !   if (iband==1) write(100+iband,*) "===================================="
1019 !   if (iband==1) write(200+iband,*) "===================================="
1020 !   if (iband==1) write(300+iband,*) "===================================="
1021 !   if (iband==1) write(400+iband,*) "===================================="
1022 !   write(100+iband,*) cwaveprj0(1,1)%cp(1,1)**2+cwaveprj0(1,1)%cp(2,1)**2 + cwaveprj0(1,2)%cp(1,1)**2+cwaveprj0(1,2)%cp(2,1)**2,&
1023 ! &                    cwaveprj0(1,1)%cp(1,2)**2+cwaveprj0(1,1)%cp(2,2)**2 + cwaveprj0(1,2)%cp(1,2)**2+cwaveprj0(1,2)%cp(2,2)**2,&
1024 ! &                    cwaveprj0(1,1)%cp(1,3)**2+cwaveprj0(1,1)%cp(2,3)**2 + cwaveprj0(1,2)%cp(1,3)**2+cwaveprj0(1,2)%cp(2,3)**2,&
1025 ! &                    cwaveprj0(1,1)%cp(1,4)**2+cwaveprj0(1,1)%cp(2,4)**2 + cwaveprj0(1,2)%cp(1,4)**2+cwaveprj0(1,2)%cp(2,4)**2,&
1026 ! &                    cwaveprj0(1,1)%cp(1,5)**2+cwaveprj0(1,1)%cp(2,5)**2 + cwaveprj0(1,2)%cp(1,5)**2+cwaveprj0(1,2)%cp(2,5)**2,&
1027 ! &                    cwaveprj0(1,1)%cp(1,6)**2+cwaveprj0(1,1)%cp(2,6)**2 + cwaveprj0(1,2)%cp(1,6)**2+cwaveprj0(1,2)%cp(2,6)**2,&
1028 ! &                    cwaveprj0(1,1)%cp(1,7)**2+cwaveprj0(1,1)%cp(2,7)**2 + cwaveprj0(1,2)%cp(1,7)**2+cwaveprj0(1,2)%cp(2,7)**2,&
1029 ! &                    cwaveprj0(1,1)%cp(1,8)**2+cwaveprj0(1,1)%cp(2,8)**2 + cwaveprj0(1,2)%cp(1,8)**2+cwaveprj0(1,2)%cp(2,8)**2
1030 !   write(200+iband,*) cwaveprj1(1,1)%cp(1,1)**2+cwaveprj1(1,1)%cp(2,1)**2 + cwaveprj1(1,2)%cp(1,1)**2+cwaveprj1(1,2)%cp(2,1)**2,&
1031 ! &                    cwaveprj1(1,1)%cp(1,2)**2+cwaveprj1(1,1)%cp(2,2)**2 + cwaveprj1(1,2)%cp(1,2)**2+cwaveprj1(1,2)%cp(2,2)**2,&
1032 ! &                    cwaveprj1(1,1)%cp(1,3)**2+cwaveprj1(1,1)%cp(2,3)**2 + cwaveprj1(1,2)%cp(1,3)**2+cwaveprj1(1,2)%cp(2,3)**2,&
1033 ! &                    cwaveprj1(1,1)%cp(1,4)**2+cwaveprj1(1,1)%cp(2,4)**2 + cwaveprj1(1,2)%cp(1,4)**2+cwaveprj1(1,2)%cp(2,4)**2,&
1034 ! &                    cwaveprj1(1,1)%cp(1,5)**2+cwaveprj1(1,1)%cp(2,5)**2 + cwaveprj1(1,2)%cp(1,5)**2+cwaveprj1(1,2)%cp(2,5)**2,&
1035 ! &                    cwaveprj1(1,1)%cp(1,6)**2+cwaveprj1(1,1)%cp(2,6)**2 + cwaveprj1(1,2)%cp(1,6)**2+cwaveprj1(1,2)%cp(2,6)**2,&
1036 ! &                    cwaveprj1(1,1)%cp(1,7)**2+cwaveprj1(1,1)%cp(2,7)**2 + cwaveprj1(1,2)%cp(1,7)**2+cwaveprj1(1,2)%cp(2,7)**2,&
1037 ! &                    cwaveprj1(1,1)%cp(1,8)**2+cwaveprj1(1,1)%cp(2,8)**2 + cwaveprj1(1,2)%cp(1,8)**2+cwaveprj1(1,2)%cp(2,8)**2
1038 !   write(300+iband,*) cwave0(1,1)**2+cwave0(2,1)**2 + cwave0(1,npw_k+1)**2+cwave0(2,npw_k+1)**2,&
1039 !                      cwave0(1,5)**2+cwave0(2,5)**2 + cwave0(1,npw_k+5)**2+cwave0(2,npw_k+5)**2,&
1040 !                      cwave0(1,9)**2+cwave0(2,9)**2 + cwave0(1,npw_k+9)**2+cwave0(2,npw_k+9)**2
1041 !   write(400+iband,*) cwave1(1,1)**2+cwave1(2,1)**2 + cwave1(1,npw_k+1)**2+cwave1(2,npw_k+1)**2,&
1042 ! &                    cwave1(1,5)**2+cwave1(2,5)**2 + cwave1(1,npw_k+5)**2+cwave1(2,npw_k+5)**2,&
1043 ! &                    cwave1(1,9)**2+cwave1(2,9)**2 + cwave1(1,npw_k+9)**2+cwave1(2,npw_k+9)**2
1044 !  end if
1045 ! end if
1046      call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj0,cwaveprj1,ipert,isppol,&
1047 &     my_natom,natom,nspinor,occ_k(iband),option_rhoij,pawrhoij1,usetimerev,wtk_k,&
1048 &     comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1049 !TEST
1050 !if (iband==1) then
1051 !  if (nspinor==1) then
1052 !   if (iband==1) v1=zero
1053 !   if (iband==1) v2=zero
1054 !   i1=1;i2=1
1055 !   v1=v1+(cwaveprj0(1,1)%cp(1,i1)*cwaveprj1(1,1)%cp(1,i2)+cwaveprj0(1,1)%cp(1,i2)*cwaveprj1(1,1)%cp(1,i1) &
1056 ! &      +cwaveprj0(1,1)%cp(2,i1)*cwaveprj1(1,1)%cp(2,i2)+cwaveprj0(1,1)%cp(2,i2)*cwaveprj1(1,1)%cp(2,i1))&
1057 ! &     *occ_k(iband)*wtk_k
1058 !   i1=1;i2=3
1059 !   v2=v2+(cwaveprj0(1,1)%cp(1,i1)*cwaveprj1(1,1)%cp(1,i2)+cwaveprj0(1,1)%cp(1,i2)*cwaveprj1(1,1)%cp(1,i1) &
1060 ! &     +cwaveprj0(1,1)%cp(2,i1)*cwaveprj1(1,1)%cp(2,i2)+cwaveprj0(1,1)%cp(2,i2)*cwaveprj1(1,1)%cp(2,i1))&
1061 ! &     *occ_k(iband)*wtk_k
1062 !   write(500,*) "A: ",v1,v2
1063 !   write(500,*) "B: ",pawrhoij1(1)%rhoij_(1,1),pawrhoij1(1)%rhoij_(4,1)
1064 !  end if
1065 !  if (nspinor==2) then
1066 !   if (iband==1) v1=zero
1067 !   if (iband==1) v2=zero
1068 !   i1=1;i2=1
1069 !   v1=v1+(cwaveprj0(1,1)%cp(1,i1)*cwaveprj1(1,1)%cp(1,i2)+cwaveprj0(1,1)%cp(1,i2)*cwaveprj1(1,1)%cp(1,i1) &
1070 ! &       +cwaveprj0(1,1)%cp(2,i1)*cwaveprj1(1,1)%cp(2,i2)+cwaveprj0(1,1)%cp(2,i2)*cwaveprj1(1,1)%cp(2,i1) &
1071 ! &       +cwaveprj0(1,2)%cp(1,i1)*cwaveprj1(1,2)%cp(1,i2)+cwaveprj0(1,2)%cp(1,i2)*cwaveprj1(1,2)%cp(1,i1) &
1072 ! &       +cwaveprj0(1,2)%cp(2,i1)*cwaveprj1(1,2)%cp(2,i2)+cwaveprj0(1,2)%cp(2,i2)*cwaveprj1(1,2)%cp(2,i1))&
1073 ! &     *occ_k(iband)*wtk_k
1074 !   i1=1;i2=3
1075 !   r1=cwaveprj0(1,1)%cp(1,i1)*cwaveprj1(1,1)%cp(1,i2)+cwaveprj0(1,1)%cp(1,i2)*cwaveprj1(1,1)%cp(1,i1) &
1076 ! &   +cwaveprj0(1,1)%cp(2,i1)*cwaveprj1(1,1)%cp(2,i2)+cwaveprj0(1,1)%cp(2,i2)*cwaveprj1(1,1)%cp(2,i1)
1077 !   r2=cwaveprj0(1,2)%cp(1,i1)*cwaveprj1(1,2)%cp(1,i2)+cwaveprj0(1,2)%cp(1,i2)*cwaveprj1(1,2)%cp(1,i1) &
1078 ! &   +cwaveprj0(1,2)%cp(2,i1)*cwaveprj1(1,2)%cp(2,i2)+cwaveprj0(1,2)%cp(2,i2)*cwaveprj1(1,2)%cp(2,i1)
1079 !  v2=v2+(cwaveprj0(1,1)%cp(1,i1)*cwaveprj1(1,1)%cp(1,i2)+cwaveprj0(1,1)%cp(1,i2)*cwaveprj1(1,1)%cp(1,i1) &
1080 ! &       +cwaveprj0(1,1)%cp(2,i1)*cwaveprj1(1,1)%cp(2,i2)+cwaveprj0(1,1)%cp(2,i2)*cwaveprj1(1,1)%cp(2,i1) &
1081 ! &       +cwaveprj0(1,2)%cp(1,i1)*cwaveprj1(1,2)%cp(1,i2)+cwaveprj0(1,2)%cp(1,i2)*cwaveprj1(1,2)%cp(1,i1) &
1082 ! &       +cwaveprj0(1,2)%cp(2,i1)*cwaveprj1(1,2)%cp(2,i2)+cwaveprj0(1,2)%cp(2,i2)*cwaveprj1(1,2)%cp(2,i1))&
1083 ! &       *occ_k(iband)*wtk_k
1084 !
1085 ! i1=1;i2=3
1086 !   write(500,*) "A: ",v2,r1,r2,&
1087 ! & cwaveprj0(1,1)%cp(1:2,i1),cwaveprj0(1,1)%cp(1:2,i2),&
1088 ! & cwaveprj0(1,2)%cp(1:2,i1),cwaveprj0(1,2)%cp(1:2,i2),&
1089 ! & cwaveprj1(1,1)%cp(1:2,i1),cwaveprj1(1,1)%cp(1:2,i2),&
1090 ! & cwaveprj1(1,2)%cp(1:2,i1),cwaveprj1(1,2)%cp(1:2,i2)
1091 !   write(500,*) "B: ",pawrhoij1(1)%rhoij_(4,1)
1092 !  end if
1093 !end if
1094    else
1095      ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,nspinor))
1096      call pawcprj_alloc(cwaveprj_tmp,ncpgr,gs_hamkq%dimcprj)
1097      choice=2
1098      call getcprj(choice,0,cwave0,cwaveprj_tmp,&
1099 &     gs_hamkq%ffnl_k,idir,gs_hamkq%indlmn,gs_hamkq%istwf_k,&
1100 &     gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,gs_hamkq%lmnmax,&
1101 &     gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
1102 &     gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,&
1103 &     gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm)
1104      call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj_tmp,cwaveprj1,ipert,isppol,&
1105 &     my_natom,gs_hamkq%natom,nspinor,occ_k(iband),option_rhoij,pawrhoij1,usetimerev,wtk_k, &
1106 &     comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1107      call pawcprj_free(cwaveprj_tmp)
1108      ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
1109    end if
1110 
1111  end if
1112 
1113  DBG_EXIT("COLL")
1114 
1115 end subroutine dfpt_accrho

ABINIT/dfpt_mkrho [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_mkrho

FUNCTION

 Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3
 from input RF and GS wavefunctions, band occupations, and k point weights.

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=wf in G space
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=first-order wf in G space
  cplex=1 if rhor1 is real, 2 if rhor1 is complex
  gprimd(3,3)=dimensional reciprocal space primitive translations
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates, GS data.
  kg1(3,mpw1*mkmem1)=reduced planewave coordinates, RF data.
  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem=Number of k points treated by this node (GS data)
  mk1mem=Number of k points treated by this node (RF data)
  mpi_enreg=information about MPI parallelization
  mpw=maximum allowed value for npw (GS wfs)
  mpw1=maximum allowed value for npw1 (RF data)
  nband_rbz(nkpt_rbz*nsppol)=number of bands to be included in summation
   at each k point for each spin channel.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
    see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt_rbz=number of k points in the reduced Brillouin zone
  npwarr(nkpt_rbz)=number of planewaves and boundary planewaves at k points
  npwar1(nkpt_rbz)=number of planewaves and boundary planewaves at k+q points
  nspden=number of spin-density components
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of symmetry elements in group (at least 1 for identity)
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation numbers for each band
   (usually 2.0) at each k point of the reduced Brillouin zone
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  rprimd(3,3)=dimensional real space primitive translations
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symrel(3,3,nsym)=symmetry matrices in real space (integers)
  ucvol=unit cell volume (Bohr**3)
  wtk_rbz(nkpt_rbz)=k point weights (they sum to 1.0).

OUTPUT

  rhog1(2,nfft)=total electron density in G space
  rhor1(cplex*nfft,nspden)=electron density in r space
   (if spin polarized, array contains total density in first half and
    spin-up density in second half)

PARENTS

      dfpt_looppert

CHILDREN

      cg_zcopy,fftpac,fourwf,sphereboundary,symrhg,timab,wrtout,xmpi_sum

SOURCE

121 subroutine dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon,istwfk_rbz,&
122 & kg,kg1,mband,mgfft,mkmem,mk1mem,mpi_enreg,mpw,mpw1,nband_rbz,&
123 & nfft,ngfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,nsppol,nsym,&
124 & occ_rbz,paral_kgb,phnons,rhog1,rhor1,rprimd,symafm,symrel,ucvol,wtk_rbz)
125 
126 
127 !This section has been created automatically by the script Abilint (TD).
128 !Do not modify the following lines by hand.
129 #undef ABI_FUNC
130 #define ABI_FUNC 'dfpt_mkrho'
131 !End of the abilint section
132 
133  implicit none
134 
135 !Arguments ------------------------------------
136 !scalars
137  integer,intent(in) :: cplex,mband,mgfft,mk1mem,mkmem,mpw,mpw1,nfft,nkpt_rbz
138  integer,intent(in) :: nspden,nspinor,nsppol,nsym,paral_kgb
139  real(dp),intent(in) :: ucvol
140  type(MPI_type),intent(in) :: mpi_enreg
141 !arrays
142  integer,intent(in) :: irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
143  integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
144  integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18),npwar1(nkpt_rbz)
145  integer,intent(in) :: npwarr(nkpt_rbz),symafm(nsym),symrel(3,3,nsym)
146  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
147  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol),gprimd(3,3)
148  real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
149  real(dp),intent(in) :: phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
150  real(dp),intent(in) :: rprimd(3,3),wtk_rbz(nkpt_rbz)
151  real(dp),intent(out) :: rhog1(2,nfft),rhor1(cplex*nfft,nspden)
152 
153 !Local variables-------------------------------
154 !scalars
155  integer,parameter :: tim_fourwf7=7,tim_rwwf15=15
156  integer,save :: nskip=0
157  integer :: bdtot_index,i1,i2,i3,iband,icg,icg1,ierr,ifft,ikg,ptr
158  integer :: ikg1,ikpt,ispden,ispinor,isppol,istwf_k,ptr1,ptr2
159  integer :: me,n1,n2,n3,n4,n5,n6,nband_k,npw1_k
160  integer :: npw_k,spaceworld
161  real(dp) :: im0,im1,re0,re1,weight
162  real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down
163  character(len=500) :: message
164 !arrays
165  integer,allocatable :: gbound(:,:),gbound1(:,:),kg1_k(:,:)
166  integer,allocatable :: kg_k(:,:)
167  real(dp) :: tsec(2)
168  real(dp),allocatable,target :: cwavef(:,:),cwavef1(:,:)
169  real(dp),allocatable :: dummy(:,:),rhoaug(:,:,:,:)
170  real(dp),allocatable :: rhoaug1(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:)
171  real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:)
172  real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:)
173  real(dp),allocatable :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:)
174 
175 ! *************************************************************************
176 
177 !DBG_ENTER("COLL")
178 
179  if(nspden==4)then
180 !  NOTE : see mkrho for the modifications needed for non-collinear treatment
181    write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
182 &   ' dfpt_mkrho : WARNING -',ch10,&
183 &   '  Linear-response calculations are under construction with nspden=4',ch10,&
184 &   ' Action : modify value of nspden in input file unless you know what you are doing.'
185 !   call wrtout(ab_out,message,'COLL')
186    call wrtout(std_out,message,'COLL')
187  end if
188 
189 !Init spaceworld
190  spaceworld=mpi_enreg%comm_cell
191  me=mpi_enreg%me_kpt
192 
193 !zero the charge density array in real space
194 !$OMP PARALLEL DO
195  do ispden=1,nspden
196    do ifft=1,cplex*nfft
197      rhor1(ifft,ispden)=zero
198    end do
199  end do
200 
201 !start loop over spin and k points
202  bdtot_index=0; icg=0; icg1=0
203 
204  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
205  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
206 
207 !Note that the dimensioning of cwavef and cwavef1 does not include nspinor
208  ABI_ALLOCATE(cwavef,(2,mpw))
209  ABI_ALLOCATE(cwavef1,(2,mpw1))
210 !Actually, rhoaug is not needed, except for strong dimensioning requirement
211  ABI_ALLOCATE(dummy,(2,1))
212  ABI_ALLOCATE(rhoaug,(n4,n5,n6,nspinor**2))
213  ABI_ALLOCATE(rhoaug1,(cplex*n4,n5,n6,nspinor**2))
214  ABI_ALLOCATE(wfraug,(2,n4,n5,n6))
215  ABI_ALLOCATE(wfraug1,(2,n4,n5,n6))
216 
217 ! EB FR Separate collinear and non-collinear magnetism
218  if (nspden /= 4) then  ! EB FR nspden check
219    do isppol=1,nsppol
220 
221      ikg=0; ikg1=0
222 
223      rhoaug1(:,:,:,:)=zero
224 
225      do ikpt=1,nkpt_rbz
226 
227        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
228        istwf_k=istwfk_rbz(ikpt)
229        npw_k=npwarr(ikpt)
230        npw1_k=npwar1(ikpt)
231 
232        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
233          bdtot_index=bdtot_index+nband_k
234          cycle
235        end if
236 
237        ABI_ALLOCATE(gbound,(2*mgfft+8,2))
238        ABI_ALLOCATE(kg_k,(3,npw_k))
239        ABI_ALLOCATE(gbound1,(2*mgfft+8,2))
240        ABI_ALLOCATE(kg1_k,(3,npw1_k))
241 
242        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
243        call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
244 
245        kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
246        call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k)
247 
248 !    Loop over bands to fft and square for rho(r)
249        do iband=1,nband_k
250          if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
251 !      Only treat occupied states
252          if (abs(occ_rbz(iband+bdtot_index))>tol8) then
253 !        Treat separately the two spinor components
254            do ispinor=1,nspinor
255 !          Obtain Fourier transform in fft box and accumulate the density
256              ptr = 1 + (ispinor-1)*npw_k + (iband-1)*npw_k*nspinor + icg
257              call cg_zcopy(npw_k, cg(1,ptr), cwavef)
258 
259 !      In these two calls, rhoaug, rhoaug1 and weight are dummy variables, and are not modified
260              call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
261 &             istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,0,paral_kgb,tim_fourwf7,weight,weight)
262 
263 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4
264              ptr = 1 + (ispinor-1)*npw1_k + (iband-1)*npw1_k*nspinor + icg1
265              call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1)
266 
267              call fourwf(cplex,rhoaug1,cwavef1,dummy,wfraug1,gbound1,gbound1,&
268 &             istwf_k,kg1_k,kg1_k,mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,&
269 &             paral_kgb,tim_fourwf7,weight,weight)
270 
271 !          Compute the weight, note that the factor 2 is
272 !          not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
273              weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol
274 
275 !          Accumulate density
276              if(cplex==2)then
277 !$OMP PARALLEL DO PRIVATE(im0,im1,re0,re1)
278                do i3=1,n3
279                  do i2=1,n2
280                    do i1=1,n1
281                      re0=wfraug(1,i1,i2,i3) ; im0=wfraug(2,i1,i2,i3)
282                      re1=wfraug1(1,i1,i2,i3); im1=wfraug1(2,i1,i2,i3)
283                      rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1)
284                      rhoaug1(2*i1  ,i2,i3,1)=rhoaug1(2*i1  ,i2,i3,1)+weight*(re0*im1-im0*re1)
285                    end do
286                  end do
287                end do
288              else
289 !$OMP PARALLEL DO
290                do i3=1,n3
291                  do i2=1,n2
292                    do i1=1,n1
293                      rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+&
294 &                     weight*( wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) + wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3)  )
295                    end do
296                  end do
297                end do
298              end if ! cplex
299            end do ! ispinor
300          else !abs(occ_rbz(iband+bdtot_index))>tol8
301            nskip=nskip+1  ! if the state is not occupied. Accumulate the number of one-way 3D ffts skipped
302          end if ! abs(occ_rbz(iband+bdtot_index))>tol8
303        end do ! iband
304 
305        ABI_DEALLOCATE(gbound)
306        ABI_DEALLOCATE(kg_k)
307        ABI_DEALLOCATE(gbound1)
308        ABI_DEALLOCATE(kg1_k)
309 
310        bdtot_index=bdtot_index+nband_k
311 
312        icg=icg+npw_k*nband_k
313        ikg=ikg+npw_k
314        icg1=icg1+npw1_k*nband_k
315        ikg1=ikg1+npw1_k
316 
317      end do ! ikpt
318 
319      if (xmpi_paral==0) then !  Write the number of one-way 3D ffts skipped until now
320        write(message,'(a,i8)')' mkrho3 : number of one-way 3D ffts skipped in mkrho3 until now =',nskip
321        call wrtout(std_out,message,'PERS')
322      end if
323 
324 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1
325 
326 !  Transfer density on augmented fft grid to normal fft grid in real space
327 !  Take also into account the spin, to place it correctly in rhor1.
328 !  Note the use of cplex
329      call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1)
330 
331    end do ! loop over isppol spins
332 
333  else ! nspden = 4
334 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335 ! Part added for the non collinear magnetism
336 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
337 ! The same lines of code are in 72_response/accrho3.F90
338 ! TODO: merge these lines in a single routine??!!
339 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
340 
341    ikg=0; ikg1=0
342 
343    rhoaug1(:,:,:,:)=zero
344 
345    do ikpt=1,nkpt_rbz
346 
347      nband_k=nband_rbz(ikpt)
348      istwf_k=istwfk_rbz(ikpt)
349      npw_k=npwarr(ikpt)
350      npw1_k=npwar1(ikpt)
351 
352      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)) then
353        bdtot_index=bdtot_index+nband_k
354        cycle
355      end if
356 
357      ABI_ALLOCATE(gbound,(2*mgfft+8,2))
358      ABI_ALLOCATE(kg_k,(3,npw_k))
359      ABI_ALLOCATE(gbound1,(2*mgfft+8,2))
360      ABI_ALLOCATE(kg1_k,(3,npw1_k))
361 
362      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
363      call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
364 
365      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
366      call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k)
367 
368 !    Loop over bands to fft and square for rho(r)
369      do iband=1,nband_k
370 
371        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,1,me)) cycle
372 
373 !      Only treat occupied states
374        if (abs(occ_rbz(iband+bdtot_index))>tol8) then
375 
376 ! Build the four components of rho. We use only norm quantities and, so fourwf.
377 
378          ABI_ALLOCATE(wfraug_up,(2,n4,n5,n6))
379          ABI_ALLOCATE(wfraug_down,(2,n4,n5,n6))
380          ABI_ALLOCATE(wfraug1_up,(2,n4,n5,n6))
381          ABI_ALLOCATE(wfraug1_down,(2,n4,n5,n6))
382          ABI_ALLOCATE(cwave0_up,(2,mpw))
383          ABI_ALLOCATE(cwave0_down,(2,mpw))
384          ABI_ALLOCATE(cwave1_up,(2,mpw1))
385          ABI_ALLOCATE(cwave1_down,(2,mpw1))
386 
387 ! EB FR build spinorial wavefunctions
388 ! Obtain Fourier transform in fft box and accumulate the density
389 ! EB FR How do we manage the following lines for non-collinear????
390 ! zero order up and down spins
391          ptr1 = 1 + (iband-1)*npw_k*nspinor + icg
392          call cg_zcopy(npw_k, cg(1,ptr1), cwave0_up)
393          ptr2 = 1 + npw_k + (iband-1)*npw_k*nspinor + icg
394          call cg_zcopy(npw_k, cg(1,ptr2), cwave0_down)
395 ! first order up and down spins
396          ptr1 = 1 + (iband-1)*npw1_k*nspinor + icg1
397          call cg_zcopy(npw1_k, cg1(1,ptr1), cwave1_up)
398          ptr2 = 1 + npw1_k + (iband-1)*npw1_k*nspinor + icg1
399          call cg_zcopy(npw1_k, cg1(1,ptr2), cwave1_down)
400 
401 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4
402 !        ptr = 1 + (ispinor-1)*npw1_k + (iband-1)*npw1_k*nspinor + icg1
403 !        call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1)
404 
405 ! EB FR lines to be managed (?????)
406 
407 ! zero order
408 !        cwave0_up => cwavef(:,1:npw_k)
409 !        cwave0_down => cwavef(:,1+npw_k:2*npw_k)
410 ! first order
411 !        cwave1_up => cwavef1(:,1:npw_k)
412 !        cwave1_down => cwavef1(:,1+npw_k:2*npw_k)
413 
414 !    The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) ?? [[cite:Gonze1997]])
415          weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol
416 !density components
417 !GS wfk Fourrier Tranform
418 ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument
419          call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gbound,gbound,istwf_k,kg_k,kg_k,&
420 &         mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
421 &         0,paral_kgb,tim_fourwf7,weight,weight)
422          call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gbound,gbound,istwf_k,kg_k,kg_k,&
423 &         mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
424 &         0,paral_kgb,tim_fourwf7,weight,weight)
425  !1st order wfk Fourrier Transform
426          call fourwf(1,rhoaug1(:,:,:,2),cwave1_up,dummy,wfraug1_up,gbound,gbound,istwf_k,kg_k,kg_k,&
427 &         mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
428 &         0,paral_kgb,tim_fourwf7,weight,weight)
429          call fourwf(1,rhoaug1(:,:,:,2),cwave1_down,dummy,wfraug1_down,gbound,gbound,istwf_k,kg_k,kg_k,&
430 &         mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
431 &         0,paral_kgb,tim_fourwf7,weight,weight)
432 
433 !    Accumulate 1st-order density (x component)
434          if (cplex==2) then
435            do i3=1,n3
436              do i2=1,n2
437                do i1=1,n1
438                  re0_up=wfraug_up(1,i1,i2,i3)  ;     im0_up=wfraug_up(2,i1,i2,i3)
439                  re1_up=wfraug1_up(1,i1,i2,i3) ;     im1_up=wfraug1_up(2,i1,i2,i3)
440                  re0_down=wfraug_down(1,i1,i2,i3)  ; im0_down=wfraug_down(2,i1,i2,i3)
441                  re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
442                  rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) !n_upup
443                  rhoaug1(2*i1  ,i2,i3,1)=zero ! imag part of n_upup at k
444                  rhoaug1(2*i1-1,i2,i3,4)=rhoaug1(2*i1-1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
445                  rhoaug1(2*i1  ,i2,i3,4)=zero ! imag part of n_dndn at k
446                  rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down &
447 &                 +im0_up*im1_down+im0_down*im1_up) ! mx; the factor two is inside weight
448                  rhoaug1(2*i1  ,i2,i3,2)=zero ! imag part of mx
449                  rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down &
450 &                 +re0_up*im1_down-im0_up*re1_down) ! my; the factor two is inside weight
451                  rhoaug1(2*i1  ,i2,i3,3)=zero ! imag part of my at k
452                end do
453              end do
454            end do
455          else
456            do i3=1,n3
457              do i2=1,n2
458                do i1=1,n1
459                  re0_up=wfraug_up(1,i1,i2,i3)  ;     im0_up=wfraug_up(2,i1,i2,i3)
460                  re1_up=wfraug1_up(1,i1,i2,i3) ;     im1_up=wfraug1_up(2,i1,i2,i3)
461                  re0_down=wfraug_down(1,i1,i2,i3)  ; im0_down=wfraug_down(2,i1,i2,i3)
462                  re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
463                  rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup
464                  rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
465                  rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down &
466 &                 +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight
467                  rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down &
468 &                 +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight
469                end do
470              end do
471            end do
472          end if
473          ABI_DEALLOCATE(wfraug_up)
474          ABI_DEALLOCATE(wfraug_down)
475          ABI_DEALLOCATE(wfraug1_up)
476          ABI_DEALLOCATE(wfraug1_down)
477          ABI_DEALLOCATE(cwave0_up)
478          ABI_DEALLOCATE(cwave0_down)
479          ABI_DEALLOCATE(cwave1_up)
480          ABI_DEALLOCATE(cwave1_down)
481 
482        end if ! occupied states
483      end do ! End loop on iband
484 
485      ABI_DEALLOCATE(gbound)
486      ABI_DEALLOCATE(kg_k)
487      ABI_DEALLOCATE(gbound1)
488      ABI_DEALLOCATE(kg1_k)
489 
490      bdtot_index=bdtot_index+nband_k
491 
492      icg=icg+npw_k*nband_k
493      ikg=ikg+npw_k
494 
495      icg1=icg1+npw1_k*nband_k
496      ikg1=ikg1+npw1_k
497 
498    end do ! End loop on ikpt
499 
500 
501    if (xmpi_paral==0) then !  Write the number of one-way 3D ffts skipped until now
502      write(message,'(a,i8)')' dfpt_mkrho : number of one-way 3D ffts skipped in mkrho3 until now =',nskip
503      call wrtout(std_out,message,'PERS')
504    end if
505 
506 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1
507 
508 !  Transfer density on augmented fft grid to normal fft grid in real space
509 !  Take also into account the spin, to place it correctly in rhor1.
510 !  Note the use of cplex
511    call fftpac(1,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1)
512  end if ! nspden /= 4
513 
514 !if (xmpi_paral==1) then
515 !call timab(63,1,tsec)
516 !call wrtout(std_out,'dfpt_mkrho: loop on k-points and spins done in parallel','COLL')
517 !call xmpi_barrier(spaceworld)
518 !call timab(63,2,tsec)
519 !end if
520 
521  ABI_DEALLOCATE(cwavef)
522  ABI_DEALLOCATE(cwavef1)
523  ABI_DEALLOCATE(dummy)
524  ABI_DEALLOCATE(rhoaug)
525  ABI_DEALLOCATE(rhoaug1)
526  ABI_DEALLOCATE(wfraug)
527  ABI_DEALLOCATE(wfraug1)
528 
529 !Recreate full rhor1 on all proc.
530  call timab(48,1,tsec)
531  call timab(71,1,tsec)
532  call xmpi_sum(rhor1,spaceworld,ierr)
533  call timab(71,2,tsec)
534  call timab(48,2,tsec)
535 
536  call symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfft,ngfft,nspden,nsppol,nsym,paral_kgb,phnons,&
537 & rhog1,rhor1,rprimd,symafm,symrel)
538 
539 !We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
540 !we also have the spin-up density, symmetrized, in rhor1(:,2).
541 
542 !DBG_EXIT("COLL")
543 
544 end subroutine dfpt_mkrho

ABINIT/m_dfpt_mkrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_mkrho

FUNCTION

 Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LSI, AR, MB, MT, SPr)
  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_dfpt_mkrho
28 
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_errors
33  use m_cgtools
34  use m_xmpi
35 
36  use m_time,            only : timab
37  use m_io_tools,        only : get_unit, iomode_from_fname
38  use m_fftcore,         only : sphereboundary
39  use m_fft,             only : fftpac, fourwf
40  use m_spacepar,        only : symrhg
41  use m_dtfil,           only : status
42  use m_hamiltonian,     only : gs_hamiltonian_type
43  use m_pawrhoij,        only : pawrhoij_type
44  use m_pawcprj,         only : pawcprj_type, pawcprj_alloc, pawcprj_free
45  use m_paw_occupancies, only : pawaccrhoij
46  use m_paral_atom,      only : get_my_atmtab
47  use m_mpinfo,          only : proc_distrb_cycle
48  use m_cgprj,           only : getcprj
49 
50  implicit none
51 
52  private