TABLE OF CONTENTS
ABINIT/dfpt_accrho [ 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
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 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 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.
SOURCE
597 subroutine dfpt_accrho(cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,& 598 & eloc0_k,gs_hamkq,iband,idir,ipert,isppol,kptopt,& 599 & mpi_enreg,natom,nband_k,ncpgr,npw_k,npw1_k,nspinor,occ_k,& 600 & option,pawrhoij1,rhoaug1,tim_fourwf,wf_corrected,& 601 & wtk_k,comm_atom,mpi_atmtab) 602 603 !Arguments ------------------------------------ 604 !scalars 605 integer,intent(in) :: cplex,iband,idir,ipert,isppol,kptopt,natom,nband_k 606 integer,intent(in) :: ncpgr,npw_k,npw1_k,nspinor,option,tim_fourwf,wf_corrected 607 integer,optional,intent(in) :: comm_atom 608 integer,optional,target,intent(in) :: mpi_atmtab(:) 609 real(dp),intent(in) :: wtk_k 610 real(dp),intent(out) :: eloc0_k 611 type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq 612 type(MPI_type),intent(in) :: mpi_enreg 613 !arrays 614 real(dp),intent(in),target :: cwave0(2,npw_k*nspinor),cwave1(2,npw1_k*nspinor),cwavef(2,npw1_k*nspinor) 615 real(dp),intent(in) :: occ_k(nband_k) 616 real(dp),intent(inout) :: rhoaug1(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc) 617 type(pawcprj_type),intent(in) :: cwaveprj0(natom,nspinor*gs_hamkq%usecprj) 618 type(pawcprj_type),intent(in) :: cwaveprj1(natom,nspinor*gs_hamkq%usepaw) 619 type(pawrhoij_type),intent(inout) :: pawrhoij1(:) 620 621 !Local variables------------------------------- 622 !scalars 623 integer,parameter :: level=14 624 integer :: choice,cplex_cprj,i1,i2,i3,ispinor,my_comm_atom,my_natom,n1,n2,n3,option_rhoij 625 logical :: my_atmtab_allocated,paral_atom 626 logical :: use_timerev,use_zeromag 627 real(dp) :: im0,im1,re0,re1,valuer,diag,offdiag,weight 628 real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down 629 !arrays 630 integer,pointer :: my_atmtab(:) 631 real(dp) :: dummy(2,1) 632 real(dp),allocatable :: rhoaug(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:) 633 real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:) 634 real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:) 635 real(dp),pointer :: cwavef_sp(:,:),cwavef_up(:,:),cwavef_down(:,:) 636 real(dp),pointer :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:) 637 real(dp),pointer :: vlocal(:,:,:,:)=>null() 638 type(pawcprj_type),allocatable :: cwaveprj_tmp(:,:) 639 640 ! ********************************************************************* 641 DBG_ENTER("COLL") 642 643 if (option/=1.and.option/=2.and.option/=3) return 644 645 !Initializations 646 ABI_MALLOC(rhoaug,(gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc)) 647 n1=gs_hamkq%ngfft(1);n2=gs_hamkq%ngfft(2);n3=gs_hamkq%ngfft(3) 648 if (option==2.or.option==3) eloc0_k=zero 649 if (option==2.or.option==3) vlocal => gs_hamkq%vlocal 650 651 !Loop on spinorial components 652 ! TODO : double loop on spinors for full rhoaug1 matrix if nspden =4 653 if (gs_hamkq%nvloc/=4) then ! see later EB FR 654 ABI_MALLOC(wfraug1,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 655 656 do ispinor=1,nspinor 657 658 ! Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy 659 ! -------------------------------------------------------------------------------------------- 660 661 ! Fourier transform of cwavef. Here, rhoaug is a dummy variable. 662 if (wf_corrected==0.or.option==2.or.option==3) then 663 if (ispinor==1) then 664 cwavef_sp => cwavef(:,1:npw1_k) 665 else 666 cwavef_sp => cwavef(:,1+npw1_k:2*npw1_k) 667 end if 668 !make an inverse FFT from cwavef_sp to wfraug1 669 call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 670 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 671 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 672 & weight,weight,gpu_option=gs_hamkq%gpu_option) 673 nullify(cwavef_sp) 674 end if 675 676 ! Compute contribution of this band to zero-order potential part of the 2nd-order total energy 677 ! NB: this is spinor diagonal 678 if (option==2.or.option==3) then 679 valuer=zero 680 do i3=1,n3 681 do i2=1,n2 682 do i1=1,n1 683 valuer=valuer+vlocal(i1,i2,i3,1)*(wfraug1(1,i1,i2,i3)**2+wfraug1(2,i1,i2,i3)**2) 684 end do 685 end do 686 end do 687 688 ! Local potential energy of this band 689 eloc0_k=eloc0_k+two*valuer/dble(gs_hamkq%nfft) 690 end if ! option 691 692 ! Part devoted to the accumulation of the 1st-order density 693 ! --------------------------------------------------------- 694 if (option==1.or.option==3) then 695 696 ! Compute 1st-order WF in real space 697 ! One needs the Fourier transform of cwave1. However, only the one of 698 ! cwavef is available. If cwavef and cwave1 differs, this Fourier 699 ! transform must be computed. In both case the result is in wfraug1. 700 if (wf_corrected==1) then 701 if (ispinor==1) then 702 cwavef_sp => cwave1(:,1:npw1_k) 703 else 704 cwavef_sp => cwave1(:,1+npw1_k:2*npw1_k) 705 end if 706 call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 707 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 708 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 709 & weight,weight,gpu_option=gs_hamkq%gpu_option) 710 nullify(cwavef_sp) 711 end if 712 713 ! Compute 0-order WF in real space 714 ! TODO: add loop over ispinor_prime here 715 ABI_MALLOC(wfraug,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 716 if (ispinor==1) then 717 cwavef_sp => cwave0(:,1:npw_k) 718 else 719 cwavef_sp => cwave0(:,1+npw_k:2*npw_k) 720 end if 721 call fourwf(1,rhoaug,cwavef_sp,dummy,wfraug,gs_hamkq%gbound_k,gs_hamkq%gbound_k,& 722 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 723 & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 724 & weight,weight,gpu_option=gs_hamkq%gpu_option) 725 nullify(cwavef_sp) 726 727 ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]]) 728 weight=two*occ_k(iband)*wtk_k/gs_hamkq%ucvol 729 ! Accumulate 1st-order density 730 if (cplex==2) then 731 do i3=1,n3 732 do i2=1,n2 733 do i1=1,n1 734 re0=wfraug(1,i1,i2,i3) ; im0=wfraug(2,i1,i2,i3) 735 re1=wfraug1(1,i1,i2,i3) ; im1=wfraug1(2,i1,i2,i3) 736 ! TODO: check which terms (ispinor ispinorp) enter a given element of rhoaug1 737 rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1) 738 rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0*im1-im0*re1) 739 end do 740 end do 741 end do 742 else 743 do i3=1,n3 744 do i2=1,n2 745 do i1=1,n1 746 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1) & 747 & +weight*(wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) & 748 & +wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3)) 749 end do 750 end do 751 end do 752 end if 753 ABI_FREE(wfraug) 754 end if ! option 755 756 end do ! Loop on spinorial components if nspden=1 or 2 757 758 ABI_FREE(wfraug1) 759 else ! nvloc = 4 760 ! The same lines of code are in 72_response/dfpt_mkrho.F90 761 ! TODO merge these lines in a single routine??!! 762 ABI_MALLOC(wfraug1_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 763 ABI_MALLOC(wfraug1_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 764 765 wfraug1_up(:,:,:,:)=zero 766 wfraug1_down(:,:,:,:)=zero 767 768 ! Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy 769 ! -------------------------------------------------------------------------------------------- 770 771 ! Fourier transform of cwavef. Here, rhoaug is a dummy variable. 772 if (wf_corrected==0.or.option==2.or.option==3) then 773 cwavef_up => cwavef(:,1:npw1_k) ! wfs up spin-polarized 774 call fourwf(cplex,rhoaug(:,:,:,1),cwavef_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 775 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 776 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 777 & weight,weight,gpu_option=gs_hamkq%gpu_option) 778 nullify(cwavef_up) 779 780 cwavef_down => cwavef(:,1+npw1_k:2*npw1_k) ! wfs down spin-polarized 781 call fourwf(cplex,rhoaug(:,:,:,1),cwavef_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 782 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 783 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 784 & weight,weight,gpu_option=gs_hamkq%gpu_option) 785 nullify(cwavef_down) 786 end if 787 if (option==2.or.option==3) then 788 valuer=zero 789 diag=zero 790 offdiag=zero 791 ! EB FR 2nd term in Eq. 91 PRB52,1096 [[cite:Gonze1995]] for non-collinear magnetism 792 do i3=1,n3 793 do i2=1,n2 794 do i1=1,n1 795 diag=vlocal(i1,i2,i3,1)*(wfraug1_up(1,i1,i2,i3)**2+wfraug1_up(2,i1,i2,i3)**2)& 796 & +vlocal(i1,i2,i3,2)*(wfraug1_down(1,i1,i2,i3)**2+wfraug1_down(2,i1,i2,i3)**2) 797 offdiag=(two*vlocal(i1,i2,i3,3)*((wfraug1_up(1,i1,i2,i3)*wfraug1_down(1,i1,i2,i3))+& 798 & (wfraug1_up(2,i1,i2,i3)*wfraug1_down(2,i1,i2,i3))))+& 799 & (two*vlocal(i1,i2,i3,4)*((-wfraug1_down(2,i1,i2,i3)*wfraug1_up(1,i1,i2,i3))+& 800 & (wfraug1_down(1,i1,i2,i3)*wfraug1_up(2,i1,i2,i3)))) 801 valuer=valuer+diag+offdiag 802 end do 803 end do 804 end do 805 ! Local potential energy of this band 806 eloc0_k=eloc0_k+two*valuer/dble(gs_hamkq%nfft) 807 end if ! option 808 809 ! Part devoted to the accumulation of the 1st-order density 810 ! --------------------------------------------------------- 811 812 ! first order 813 ! 814 if (option==1.or.option==3) then 815 816 !SPr: condition on wf_corrected not to do FFTs of the same Bloch functions 817 if (wf_corrected==1) then 818 cwave1_up => cwave1(:,1:npw1_k) 819 cwave1_down => cwave1(:,1+npw1_k:2*npw1_k) 820 wfraug1_up(:,:,:,:)=zero 821 wfraug1_down(:,:,:,:)=zero 822 823 call fourwf(cplex,rhoaug(:,:,:,1),cwave1_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 824 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 825 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 826 & weight,weight,gpu_option=gs_hamkq%gpu_option) 827 nullify(cwave1_up) 828 829 call fourwf(cplex,rhoaug(:,:,:,1),cwave1_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,& 830 & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 831 & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 832 & weight,weight,gpu_option=gs_hamkq%gpu_option) 833 nullify(cwave1_down) 834 end if 835 836 837 ! EB FR build spinorial wavefunctions 838 ! zero order 839 cwave0_up => cwave0(:,1:npw_k) 840 cwave0_down => cwave0(:,1+npw_k:2*npw_k) 841 ABI_MALLOC(wfraug_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 842 ABI_MALLOC(wfraug_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 843 wfraug_up(:,:,:,:)=zero 844 wfraug_down(:,:,:,:)=zero 845 ! 846 !density components 847 !GS wfk Fourrier Tranform 848 ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument 849 call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gs_hamkq%gbound_k,gs_hamkq%gbound_k,& 850 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 851 & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 852 & weight,weight,gpu_option=gs_hamkq%gpu_option) 853 nullify(cwave0_up) 854 call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gs_hamkq%gbound_k,gs_hamkq%gbound_k,& 855 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 856 & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,& 857 & weight,weight,gpu_option=gs_hamkq%gpu_option) 858 nullify(cwave0_down) 859 ! Accumulate 1st-order density (x component) 860 re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero 861 re1_down=zero;im1_down=zero 862 ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]]) 863 ! SPr: the following treatment with factor=2 is ok for perturbations not breaking the 864 ! time reversal symmetry of the Hamiltonian (due to Kramer's degeneracy) hence 865 ! not applicable for magnetic field perturbation (for phonons with SOC, H^(0) has 866 ! time reversal symmetry though). The formulas below are rectified in dfpt_scfcv 867 ! in case of broken time-reversal upon reconstructing rhor1_pq and rhor1_mq. 868 weight=two*occ_k(iband)*wtk_k/gs_hamkq%ucvol 869 if (cplex==2) then 870 do i3=1,n3 871 do i2=1,n2 872 do i1=1,n1 873 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 874 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 875 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 876 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 877 !SPr: in case of +q/-q calculation, the factor will be corrected later from dfpt_scfcv level 878 ! along with the reconstruction of correct rhor1_{+q} and rhor1_{-q} 879 ! here, rhoaug1_{sigma,sigma'} = \sum_{n,k} u1_{sigma} u0*_{sigma'} independent of the sign of q 880 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 881 rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0_up*im1_up-im0_up*re1_up) 882 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 883 rhoaug1(2*i1 ,i2,i3,4)=rhoaug1(2*i1 ,i2,i3,4)+weight*(re0_down*im1_down-im0_down*re1_down) 884 885 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] 886 & +weight*(re1_down*re0_up+im1_down*im0_up) 887 rhoaug1(2*i1 ,i2,i3,2)=rhoaug1(2*i1 ,i2,i3,2)+weight*(-re1_up*im0_down+im1_up*re0_down)& !Im[m1x] 888 & +weight*(-re1_down*im0_up+im1_down*re0_up) 889 890 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] 891 & +weight*(-re1_down*im0_up+im1_down*re0_up) 892 rhoaug1(2*i1 ,i2,i3,3)=rhoaug1(2*i1 ,i2,i3,3)+weight*(+re1_up*re0_down+im1_up*im0_down)& !Im[m1y] 893 & +weight*(-re1_down*re0_up-im1_down*im0_up) 894 end do 895 end do 896 end do 897 else !cplex 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 do i3=1,n3 901 do i2=1,n2 902 do i1=1,n1 903 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 904 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 905 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 906 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 907 908 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup 909 rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn 910 rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down & 911 & +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight 912 rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down & 913 & +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight 914 end do 915 end do 916 end do 917 end if !cplex 918 919 ABI_FREE(wfraug_up) 920 ABI_FREE(wfraug_down) 921 922 end if ! option 923 924 ABI_FREE(wfraug1_up) 925 ABI_FREE(wfraug1_down) 926 927 end if ! nvloc /= 4 928 929 ABI_FREE(rhoaug) 930 931 !Part devoted to the accumulation of the 1st-order occupation matrix in PAW case 932 ! TODO: parse for more nspden 4 dependencies on spinors 933 ! EB FR CHECK: to be modified for non-collinear????? 934 !------------------------------------------------------------------------------- 935 936 if ((option==1.or.option==3).and.gs_hamkq%usepaw==1) then 937 938 ! Set up parallelism over atoms 939 my_natom=natom; if(gs_hamkq%usepaw==1) my_natom=size(pawrhoij1) 940 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 941 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 942 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 943 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 944 945 cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1 946 option_rhoij=2 947 use_timerev=(kptopt>0.and.kptopt<3) 948 use_zeromag=.false.;if (my_natom>0) use_zeromag=(pawrhoij1(1)%nspden==4.and.gs_hamkq%nvloc==1) 949 950 if (gs_hamkq%usecprj==1) then 951 call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj0,cwaveprj1,ipert,isppol,my_natom,& 952 & natom,nspinor,occ_k(iband),option_rhoij,pawrhoij1,use_timerev,use_zeromag,wtk_k,& 953 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 954 else 955 ABI_MALLOC(cwaveprj_tmp,(natom,nspinor)) 956 call pawcprj_alloc(cwaveprj_tmp,ncpgr,gs_hamkq%dimcprj) 957 choice=2 958 call getcprj(choice,0,cwave0,cwaveprj_tmp,& 959 & gs_hamkq%ffnl_k,idir,gs_hamkq%indlmn,gs_hamkq%istwf_k,& 960 & gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,gs_hamkq%lmnmax,& 961 & gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,& 962 & gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,& 963 & gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm) 964 call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj_tmp,cwaveprj1,ipert,isppol,my_natom,& 965 & gs_hamkq%natom,nspinor,occ_k(iband),option_rhoij,pawrhoij1,use_timerev,use_zeromag,wtk_k, & 966 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 967 call pawcprj_free(cwaveprj_tmp) 968 ABI_FREE(cwaveprj_tmp) 969 end if 970 971 end if 972 973 DBG_EXIT("COLL") 974 975 end subroutine dfpt_accrho
ABINIT/dfpt_mkrho [ 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)
SOURCE
110 subroutine dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon,istwfk_rbz,& 111 & kg,kg1,mband,mband_mem,mgfft,mkmem,mk1mem,mpi_enreg,mpw,mpw1,nband_rbz,& 112 & nfft,ngfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,nsppol,nsym,& 113 & occ_rbz,phnons,rhog1,rhor1,rprimd,symafm,symrel,tnons,ucvol,wtk_rbz) 114 115 !Arguments ------------------------------------ 116 !scalars 117 integer,intent(in) :: cplex,mband,mband_mem,mgfft,mk1mem,mkmem,mpw,mpw1,nfft,nkpt_rbz 118 integer,intent(in) :: nspden,nspinor,nsppol,nsym 119 real(dp),intent(in) :: ucvol 120 type(MPI_type),intent(in) :: mpi_enreg 121 !arrays 122 integer,intent(in) :: irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)) 123 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 124 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18),npwar1(nkpt_rbz) 125 integer,intent(in) :: npwarr(nkpt_rbz),symafm(nsym),symrel(3,3,nsym) 126 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol) 127 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol),gprimd(3,3) 128 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 129 real(dp),intent(in) :: phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)) 130 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym) 131 real(dp),intent(in) :: wtk_rbz(nkpt_rbz) 132 real(dp),intent(out) :: rhog1(2,nfft),rhor1(cplex*nfft,nspden) 133 134 !Local variables------------------------------- 135 !scalars 136 integer,parameter :: tim_fourwf7=7,tim_rwwf15=15 137 integer,save :: nskip=0 138 integer :: bdtot_index,i1,i2,i3,iband,icg,icg1,ierr,ifft,ikg,ptr 139 integer :: iband_me 140 integer :: ikg1,ikpt,ispden,ispinor,isppol,istwf_k,ptr1,ptr2 141 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,npw1_k 142 integer :: npw_k,spaceworld 143 real(dp) :: im0,im1,re0,re1,weight 144 real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down 145 character(len=500) :: message 146 !arrays 147 integer,allocatable :: gbound(:,:),gbound1(:,:),kg1_k(:,:) 148 integer,allocatable :: kg_k(:,:) 149 real(dp) :: tsec(2) 150 real(dp),allocatable,target :: cwavef(:,:),cwavef1(:,:) 151 real(dp),allocatable :: dummy(:,:),rhoaug(:,:,:,:) 152 real(dp),allocatable :: rhoaug1(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:) 153 real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:) 154 real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:) 155 real(dp),allocatable :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:) 156 157 ! ************************************************************************* 158 159 !DBG_ENTER("COLL") 160 161 if(nspden==4)then 162 ! NOTE: see mkrho for the modifications needed for non-collinear treatment 163 write(message, '(3a)' )& 164 ' Linear-response calculations are under construction with nspden=4',ch10,& 165 ' Action: modify value of nspden in input file unless you know what you are doing.' 166 ABI_WARNING(message) 167 end if 168 169 !Init spaceworld 170 spaceworld=mpi_enreg%comm_cell 171 me=mpi_enreg%me_kpt 172 173 !zero the charge density array in real space 174 !$OMP PARALLEL DO 175 do ispden=1,nspden 176 do ifft=1,cplex*nfft 177 rhor1(ifft,ispden)=zero 178 end do 179 end do 180 181 !start loop over spin and k points 182 bdtot_index=0; icg=0; icg1=0 183 184 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 185 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing 186 187 !Note that the dimensioning of cwavef and cwavef1 does not include nspinor 188 ABI_MALLOC(cwavef,(2,mpw)) 189 ABI_MALLOC(cwavef1,(2,mpw1)) 190 !Actually, rhoaug is not needed, except for strong dimensioning requirement 191 ABI_MALLOC(dummy,(2,1)) 192 ABI_MALLOC(rhoaug,(n4,n5,n6,nspinor**2)) 193 ABI_MALLOC(rhoaug1,(cplex*n4,n5,n6,nspinor**2)) 194 ABI_MALLOC(wfraug,(2,n4,n5,n6)) 195 ABI_MALLOC(wfraug1,(2,n4,n5,n6)) 196 197 ! EB FR Separate collinear and non-collinear magnetism 198 if (nspden /= 4) then ! EB FR nspden check 199 do isppol=1,nsppol 200 201 ikg=0; ikg1=0 202 203 rhoaug1(:,:,:,:)=zero 204 205 do ikpt=1,nkpt_rbz 206 207 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 208 istwf_k=istwfk_rbz(ikpt) 209 npw_k=npwarr(ikpt) 210 npw1_k=npwar1(ikpt) 211 212 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 213 bdtot_index=bdtot_index+nband_k 214 cycle 215 end if 216 217 ABI_MALLOC(gbound,(2*mgfft+8,2)) 218 ABI_MALLOC(kg_k,(3,npw_k)) 219 ABI_MALLOC(gbound1,(2*mgfft+8,2)) 220 ABI_MALLOC(kg1_k,(3,npw1_k)) 221 222 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 223 call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k) 224 225 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 226 call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k) 227 228 ! Loop over bands to fft and square for rho(r) 229 iband_me = 0 230 do iband=1,nband_k 231 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 232 iband_me = iband_me + 1 233 ! Only treat occupied states 234 if (abs(occ_rbz(iband+bdtot_index))>tol8) then 235 ! Treat separately the two spinor components 236 do ispinor=1,nspinor 237 ! Obtain Fourier transform in fft box and accumulate the density 238 ptr = 1 + (ispinor-1)*npw_k + (iband_me-1)*npw_k*nspinor + icg 239 call cg_zcopy(npw_k, cg(1,ptr), cwavef) 240 241 ! In these two calls, rhoaug, rhoaug1 and weight are dummy variables, and are not modified 242 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 243 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,0,tim_fourwf7,weight,weight) 244 245 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4 246 ptr = 1 + (ispinor-1)*npw1_k + (iband_me-1)*npw1_k*nspinor + icg1 247 call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1) 248 249 call fourwf(cplex,rhoaug1,cwavef1,dummy,wfraug1,gbound1,gbound1,& 250 & istwf_k,kg1_k,kg1_k,mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,& 251 & tim_fourwf7,weight,weight) 252 253 ! Compute the weight, note that the factor 2 is 254 ! not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]]) 255 weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol 256 257 ! Accumulate density 258 if(cplex==2)then 259 !$OMP PARALLEL DO PRIVATE(im0,im1,re0,re1) 260 do i3=1,n3 261 do i2=1,n2 262 do i1=1,n1 263 re0=wfraug(1,i1,i2,i3) ; im0=wfraug(2,i1,i2,i3) 264 re1=wfraug1(1,i1,i2,i3); im1=wfraug1(2,i1,i2,i3) 265 rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1) 266 rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0*im1-im0*re1) 267 end do 268 end do 269 end do 270 else 271 !$OMP PARALLEL DO 272 do i3=1,n3 273 do i2=1,n2 274 do i1=1,n1 275 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+& 276 & weight*( wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) + wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3) ) 277 end do 278 end do 279 end do 280 end if ! cplex 281 end do ! ispinor 282 else !abs(occ_rbz(iband+bdtot_index))>tol8 283 nskip=nskip+1 ! if the state is not occupied. Accumulate the number of one-way 3D ffts skipped 284 end if ! abs(occ_rbz(iband+bdtot_index))>tol8 285 286 end do ! iband 287 288 ABI_FREE(gbound) 289 ABI_FREE(kg_k) 290 ABI_FREE(gbound1) 291 ABI_FREE(kg1_k) 292 293 bdtot_index=bdtot_index+nband_k 294 295 ! only increase indices for my bands on my proc 296 icg=icg+npw_k*mband_mem*nspinor 297 ikg=ikg+npw_k 298 299 icg1=icg1+npw1_k*mband_mem*nspinor 300 ikg1=ikg1+npw1_k 301 302 end do ! ikpt 303 304 if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now 305 write(message,'(a,i8)')' mkrho3 : number of one-way 3D ffts skipped in mkrho3 until now =',nskip 306 call wrtout(std_out,message,'PERS') 307 end if 308 309 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1 310 311 ! Transfer density on augmented fft grid to normal fft grid in real space 312 ! Take also into account the spin, to place it correctly in rhor1. 313 ! Note the use of cplex 314 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1) 315 316 end do ! loop over isppol spins 317 318 else ! nspden = 4 319 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 320 ! Part added for the non collinear magnetism 321 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 322 ! The same lines of code are in 72_response/accrho3.F90 323 ! TODO: merge these lines in a single routine??!! 324 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 325 326 ikg=0; ikg1=0 327 328 rhoaug1(:,:,:,:)=zero 329 330 do ikpt=1,nkpt_rbz 331 332 nband_k=nband_rbz(ikpt) 333 istwf_k=istwfk_rbz(ikpt) 334 npw_k=npwarr(ikpt) 335 npw1_k=npwar1(ikpt) 336 337 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)) then 338 bdtot_index=bdtot_index+nband_k 339 cycle 340 end if 341 342 ABI_MALLOC(gbound,(2*mgfft+8,2)) 343 ABI_MALLOC(kg_k,(3,npw_k)) 344 ABI_MALLOC(gbound1,(2*mgfft+8,2)) 345 ABI_MALLOC(kg1_k,(3,npw1_k)) 346 347 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 348 call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k) 349 350 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 351 call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k) 352 353 ! Loop over bands to fft and square for rho(r) 354 iband_me = 0 355 do iband=1,nband_k 356 357 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,1,me)) cycle 358 iband_me = iband_me + 1 359 360 ! Only treat occupied states 361 if (abs(occ_rbz(iband+bdtot_index))>tol8) then 362 363 ! Build the four components of rho. We use only norm quantities and, so fourwf. 364 365 ABI_MALLOC(wfraug_up,(2,n4,n5,n6)) 366 ABI_MALLOC(wfraug_down,(2,n4,n5,n6)) 367 ABI_MALLOC(wfraug1_up,(2,n4,n5,n6)) 368 ABI_MALLOC(wfraug1_down,(2,n4,n5,n6)) 369 ABI_MALLOC(cwave0_up,(2,mpw)) 370 ABI_MALLOC(cwave0_down,(2,mpw)) 371 ABI_MALLOC(cwave1_up,(2,mpw1)) 372 ABI_MALLOC(cwave1_down,(2,mpw1)) 373 374 ! EB FR build spinorial wavefunctions 375 ! Obtain Fourier transform in fft box and accumulate the density 376 ! EB FR How do we manage the following lines for non-collinear???? 377 ! zero order up and down spins 378 ptr1 = 1 + (iband_me-1)*npw_k*nspinor + icg 379 call cg_zcopy(npw_k, cg(1,ptr1), cwave0_up) 380 ptr2 = 1 + npw_k + (iband_me-1)*npw_k*nspinor + icg 381 call cg_zcopy(npw_k, cg(1,ptr2), cwave0_down) 382 ! first order up and down spins 383 ptr1 = 1 + (iband_me-1)*npw1_k*nspinor + icg1 384 call cg_zcopy(npw1_k, cg1(1,ptr1), cwave1_up) 385 ptr2 = 1 + npw1_k + (iband_me-1)*npw1_k*nspinor + icg1 386 call cg_zcopy(npw1_k, cg1(1,ptr2), cwave1_down) 387 388 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4 389 ! ptr = 1 + (ispinor-1)*npw1_k + (iband_me-1)*npw1_k*nspinor + icg1 390 ! call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1) 391 392 ! EB FR lines to be managed (?????) 393 394 ! zero order 395 ! cwave0_up => cwavef(:,1:npw_k) 396 ! cwave0_down => cwavef(:,1+npw_k:2*npw_k) 397 ! first order 398 ! cwave1_up => cwavef1(:,1:npw_k) 399 ! cwave1_down => cwavef1(:,1+npw_k:2*npw_k) 400 401 ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) ?? [[cite:Gonze1997]]) 402 weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol 403 !density components 404 !GS wfk Fourrier Tranform 405 ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument 406 call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gbound,gbound,istwf_k,kg_k,kg_k,& 407 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 408 & 0,tim_fourwf7,weight,weight) 409 call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gbound,gbound,istwf_k,kg_k,kg_k,& 410 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 411 & 0,tim_fourwf7,weight,weight) 412 !1st order wfk Fourrier Transform 413 call fourwf(1,rhoaug1(:,:,:,2),cwave1_up,dummy,wfraug1_up,gbound,gbound,istwf_k,kg_k,kg_k,& 414 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 415 & 0,tim_fourwf7,weight,weight) 416 call fourwf(1,rhoaug1(:,:,:,2),cwave1_down,dummy,wfraug1_down,gbound,gbound,istwf_k,kg_k,kg_k,& 417 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 418 & 0,tim_fourwf7,weight,weight) 419 420 ! Accumulate 1st-order density (x component) 421 if (cplex==2) then 422 do i3=1,n3 423 do i2=1,n2 424 do i1=1,n1 425 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 426 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 427 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 428 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 429 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 430 rhoaug1(2*i1 ,i2,i3,1)=zero ! imag part of n_upup at k 431 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 432 rhoaug1(2*i1 ,i2,i3,4)=zero ! imag part of n_dndn at k 433 rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down & 434 & +im0_up*im1_down+im0_down*im1_up) ! mx; the factor two is inside weight 435 rhoaug1(2*i1 ,i2,i3,2)=zero ! imag part of mx 436 rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down & 437 & +re0_up*im1_down-im0_up*re1_down) ! my; the factor two is inside weight 438 rhoaug1(2*i1 ,i2,i3,3)=zero ! imag part of my at k 439 end do 440 end do 441 end do 442 else 443 do i3=1,n3 444 do i2=1,n2 445 do i1=1,n1 446 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 447 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 448 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 449 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 450 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup 451 rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn 452 rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down & 453 & +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight 454 rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down & 455 & +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight 456 end do 457 end do 458 end do 459 end if 460 ABI_FREE(wfraug_up) 461 ABI_FREE(wfraug_down) 462 ABI_FREE(wfraug1_up) 463 ABI_FREE(wfraug1_down) 464 ABI_FREE(cwave0_up) 465 ABI_FREE(cwave0_down) 466 ABI_FREE(cwave1_up) 467 ABI_FREE(cwave1_down) 468 469 end if ! occupied states 470 end do ! End loop on iband 471 472 ABI_FREE(gbound) 473 ABI_FREE(kg_k) 474 ABI_FREE(gbound1) 475 ABI_FREE(kg1_k) 476 477 bdtot_index=bdtot_index+nband_k 478 479 ! only increase indices for my bands on my proc 480 icg=icg+npw_k*mband_mem*nspinor 481 ikg=ikg+npw_k 482 483 icg1=icg1+npw1_k*mband_mem*nspinor 484 ikg1=ikg1+npw1_k 485 486 ikg=ikg+npw_k 487 ikg1=ikg1+npw1_k 488 489 end do ! End loop on ikpt 490 491 492 if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now 493 write(message,'(a,i8)')' dfpt_mkrho : number of one-way 3D ffts skipped in mkrho3 until now =',nskip 494 call wrtout(std_out,message,'PERS') 495 end if 496 497 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1 498 499 ! Transfer density on augmented fft grid to normal fft grid in real space 500 ! Take also into account the spin, to place it correctly in rhor1. 501 ! Note the use of cplex 502 call fftpac(1,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1) 503 end if ! nspden /= 4 504 505 !if (xmpi_paral==1) then 506 !call timab(63,1,tsec) 507 !call wrtout(std_out,'dfpt_mkrho: loop on k-points and spins done in parallel','COLL') 508 !call xmpi_barrier(spaceworld) 509 !call timab(63,2,tsec) 510 !end if 511 512 ABI_FREE(cwavef) 513 ABI_FREE(cwavef1) 514 ABI_FREE(dummy) 515 ABI_FREE(rhoaug) 516 ABI_FREE(rhoaug1) 517 ABI_FREE(wfraug) 518 ABI_FREE(wfraug1) 519 520 !Recreate full rhor1 on all proc. 521 !TODO : check this sums correctly on bands as well as k 522 call timab(48,1,tsec) 523 call timab(71,1,tsec) 524 call xmpi_sum(rhor1,spaceworld,ierr) 525 call timab(71,2,tsec) 526 call timab(48,2,tsec) 527 528 call symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfft,ngfft,nspden,nsppol,nsym,phnons,& 529 rhog1,rhor1,rprimd,symafm,symrel,tnons) 530 531 !We now have both rho(r) and rho(G), symmetrized, and if nsppol=2 532 !we also have the spin-up density, symmetrized, in rhor1(:,2). 533 534 !DBG_EXIT("COLL") 535 536 end subroutine dfpt_mkrho
ABINIT/m_dfpt_mkrho [ Modules ]
NAME
m_dfpt_mkrho
FUNCTION
Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_dfpt_mkrho 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_cgtools 28 use m_xmpi 29 30 31 use defs_abitypes, only : MPI_type 32 use m_time, only : timab 33 use m_io_tools, only : get_unit, iomode_from_fname 34 use m_fftcore, only : sphereboundary 35 use m_fft, only : fftpac, fourwf 36 use m_spacepar, only : symrhg 37 use m_hamiltonian, only : gs_hamiltonian_type 38 use m_pawrhoij, only : pawrhoij_type 39 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 40 use m_paw_occupancies, only : pawaccrhoij 41 use m_paral_atom, only : get_my_atmtab 42 use m_mpinfo, only : proc_distrb_cycle 43 use m_cgprj, only : getcprj 44 45 implicit none 46 47 private