TABLE OF CONTENTS
ABINIT/fock2ACE [ Functions ]
NAME
fock2ACE
FUNCTION
Compute nonlocal contribution to the Fock part of the hamiltonian in the ACE formalism. optionally contribution to Fock forces
INPUTS
cg(2,mcg)=wavefunctions (may be read from disk file) cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> fock <type(fock_type)>= quantities to calculate Fock exact exchange istwfk(nkpt)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis kpt(3,nkpt)=k points in reduced coordinates mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mkmem=number of k points treated by this node. mpi_enreg=information about MPI parallelization mpsang= mpw= maximum number of plane waves my_natom=number of atoms treated by current processor natom=number of atoms in cell. nband(nkpt)=number of bands at each k point nfft=number of FFT grid points ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt=number of k points in Brillouin zone nloalg(3)=governs the choice of the algorithm for non-local operator. npwarr(nkpt)=number of planewaves in basis and boundary at each k nspden=Number of spin Density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points optfor=1 if computation of forces is required paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(natom)=type of each atom usecprj=1 if cprj datastructure has been allocated gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU) wtk(nkpt)=weight associated with each k point xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
OUTPUT
fock%fockACE(ikpt,isppol)%xi if optfor=1, fock%fock_common%forces
SOURCE
804 subroutine fock2ACE(cg,cprj,fock,istwfk,kg,kpt,mband,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,& 805 & mpw,my_natom,natom,nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,& 806 & ntypat,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,typat,usecprj,gpu_option,wtk,xred,ylm) 807 808 !Arguments ------------------------------------ 809 !scalars 810 integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nkpt 811 integer,intent(in) :: nspden,nsppol,nspinor,ntypat,optfor 812 integer,intent(in) :: usecprj,gpu_option 813 type(MPI_type),intent(inout) :: mpi_enreg 814 type(pseudopotential_type),intent(in) :: psps 815 !arrays 816 integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol) 817 integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt) 818 integer,intent(in) :: typat(natom) 819 real(dp),intent(in) :: cg(2,mcg) 820 real(dp),intent(in) :: kpt(3,nkpt) 821 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom) 822 real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom) 823 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 824 type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj) 825 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 826 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 827 type(fock_type),pointer, intent(inout) :: fock 828 !Local variables------------------------------- 829 !scalars 830 integer :: bandpp,bdtot_index,dimffnl,iband,iband_cprj,iband_last,ibg,icg,ider 831 integer :: ierr,info,idir,ikg,ikpt,ilm,ipw,isppol,istwf_k,kk,ll 832 integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,ndat,nkpg 833 integer :: npw_k,spaceComm 834 integer :: use_ACE_old 835 integer :: blocksize,iblock,jblock,iblocksize,jblocksize,nblockbd 836 !integer, save :: counter=0 837 type(gs_hamiltonian_type) :: gs_hamk 838 logical :: compute_gbound 839 character(len=500) :: msg 840 type(fock_common_type),pointer :: fockcommon 841 !arrays 842 integer,allocatable :: kg_k(:,:) 843 real(dp) :: kpoint(3),rmet(3,3),tsec(2) 844 real(dp),allocatable :: bb(:,:,:),cwavef(:,:),cwavefk(:,:),ffnl_sav(:,:,:,:) 845 real(dp),allocatable :: kpg_k(:,:),kpg_k_sav(:,:) 846 real(dp),allocatable :: mkl(:,:,:),occblock(:),ph3d(:,:,:),ph3d_sav(:,:,:) 847 real(dp),allocatable :: wi(:,:,:),weight(:),ylm_k(:,:),ylmgr_k(:,:,:) 848 real(dp),allocatable,target :: ffnl(:,:,:,:) 849 type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null() 850 type(pawcprj_type),target,allocatable :: cwaveprj(:,:) 851 type(pawcprj_type),pointer :: cwaveprj_idat(:,:) 852 853 !************************************************************************* 854 855 call timab(1560,1,tsec) 856 call timab(1561,1,tsec) 857 858 !DEBUG 859 !if(counter>0)return 860 !counter=counter+1 861 !ENDDEBUG 862 863 !Init mpicomm and me 864 if(mpi_enreg%paral_kgb==1)then 865 spaceComm=mpi_enreg%comm_kpt 866 me_distrb=mpi_enreg%me_kpt 867 else 868 !* In case of HF calculation 869 if (mpi_enreg%paral_hf==1) then 870 spaceComm=mpi_enreg%comm_kpt 871 me_distrb=mpi_enreg%me_kpt 872 else 873 spaceComm=mpi_enreg%comm_cell 874 me_distrb=mpi_enreg%me_cell 875 end if 876 end if 877 878 !Some initializations 879 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 880 compute_gbound=.true. 881 fockcommon => fock%fock_common 882 use_ACE_old=fockcommon%use_ACE 883 fockcommon%use_ACE=0 884 fockcommon%e_fock0=zero 885 886 !Initialize Hamiltonian (k- and spin-independent terms) 887 888 call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,& 889 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj,& 890 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 891 & paw_ij=paw_ij,ph1d=ph1d,fock=fock,& 892 & gpu_option=gpu_option) 893 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 894 fockcommon%use_ACE=use_ACE_old 895 896 !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted) 897 if (psps%usepaw==1) then 898 call pawcprj_reorder(cprj,gs_hamk%atindx) 899 end if 900 901 !LOOP OVER SPINS 902 bdtot_index=0;ibg=0;icg=0 903 call timab(1561,2,tsec) ; call timab(1562,-1,tsec) 904 905 do isppol=1,nsppol 906 fockcommon%isppol=isppol 907 ! Continue to initialize the Hamiltonian (PAW DIJ coefficients) 908 call gs_hamk%load_spin(isppol,with_nonlocal=.true.) 909 910 ! Loop over k points 911 ikg=0 912 do ikpt=1,nkpt 913 fockcommon%ikpt=ikpt 914 nband_k=nband(ikpt+(isppol-1)*nkpt) 915 npw_k=npwarr(ikpt) 916 kpoint(:)=kpt(:,ikpt) 917 istwf_k=istwfk(ikpt) 918 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 919 bdtot_index=bdtot_index+nband_k 920 cycle 921 end if 922 923 ! Parallelism over FFT and/or bands: define sizes and tabs 924 if (mpi_enreg%paral_kgb==1) then 925 my_ikpt=mpi_enreg%my_kpttab(ikpt) 926 nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp) 927 bandpp=mpi_enreg%bandpp 928 my_bandfft_kpt => bandfft_kpt(my_ikpt) 929 else 930 my_ikpt=ikpt 931 bandpp=mpi_enreg%bandpp 932 nblockbd=nband_k/bandpp 933 end if 934 blocksize=nband_k/nblockbd 935 mband_cprj=mband/mpi_enreg%nproc_band 936 nband_cprj_k=nband_k/mpi_enreg%nproc_band 937 938 ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize)) 939 if (psps%usepaw==1) then 940 ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp)) 941 call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj) 942 else 943 ABI_MALLOC(cwaveprj,(0,0)) 944 end if 945 946 ABI_MALLOC(kg_k,(3,mpw)) 947 !$OMP PARALLEL DO 948 do ipw=1,npw_k 949 kg_k(:,ipw)=kg(:,ipw+ikg) 950 end do 951 952 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 953 ABI_MALLOC(ylmgr_k,(0,0,0)) 954 if (psps%useylm==1) then 955 !$OMP PARALLEL DO COLLAPSE(2) 956 do ilm=1,mpsang*mpsang 957 do ipw=1,npw_k 958 ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm) 959 end do 960 end do 961 end if 962 963 ! Compute (k+G) vectors 964 nkpg=3*nloalg(3) 965 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 966 if (nkpg>0) then 967 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 968 end if 969 970 971 ! Compute nonlocal form factors ffnl at all (k+G) 972 ider=0;idir=0;dimffnl=1 973 ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 974 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 975 & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,& 976 & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 977 978 ! Load k-dependent part in the Hamiltonian datastructure 979 ! - Compute 3D phase factors 980 ! - Prepare various tabs in case of band-FFT parallelism 981 ! - Load k-dependent quantities in the Hamiltonian 982 983 ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk)) 984 call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,& 985 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.) 986 987 ! Load band-FFT tabs (transposed k-dependent arrays) 988 if (mpi_enreg%paral_kgb==1) then 989 call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav) 990 call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg) 991 call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, & 992 & kg_k =my_bandfft_kpt%kg_k_gather, & 993 & kpg_k =my_bandfft_kpt%kpg_k_gather, & 994 ffnl_k =my_bandfft_kpt%ffnl_gather, & 995 ph3d_k =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound) 996 end if 997 998 ! The following is now wrong. In sequential, nblockbd=nband_k/bandpp 999 ! blocksize= bandpp (JB 2016/04/16) 1000 ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1 1001 ! 1002 ABI_MALLOC(occblock,(blocksize)) 1003 ABI_MALLOC(weight,(blocksize)) 1004 occblock=zero;weight=zero 1005 1006 if (fockcommon%optfor) then 1007 fockcommon%forces_ikpt=zero 1008 end if 1009 1010 ABI_MALLOC(wi,(2,npw_k*my_nspinor*blocksize,nblockbd)) 1011 wi=zero 1012 ABI_MALLOC(mkl,(2,nband_k,nband_k)) 1013 mkl=zero 1014 ! Calculate all the Wi for the current k-point 1015 1016 do iblock=1,nblockbd 1017 1018 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k) 1019 iband_cprj=(iblock-1)*bandpp+1 1020 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle 1021 1022 ! Select occupied bandsddk 1023 occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 1024 weight(:)=wtk(ikpt)*occblock(:) 1025 1026 ! Load contribution from n,k 1027 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 1028 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 1029 if (psps%usepaw==1) then 1030 call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,& 1031 & mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,& 1032 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 1033 end if 1034 1035 if (mpi_enreg%paral_kgb==1) then 1036 msg='fock2ACE: Paral_kgb is not yet implemented for fock calculations' 1037 ABI_BUG(msg) 1038 end if 1039 ndat=mpi_enreg%bandpp 1040 if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj 1041 1042 do iblocksize=1,blocksize 1043 fockcommon%ieigen=(iblock-1)*blocksize+iblocksize 1044 fockcommon%iband=(iblock-1)*blocksize+iblocksize 1045 if (gs_hamk%usepaw==1) then 1046 cwaveprj_idat => cwaveprj(:,(iblocksize-1)*my_nspinor+1:iblocksize*my_nspinor) 1047 end if 1048 call timab(1562,2,tsec) ; call timab(1563,-1,tsec) 1049 call fock_getghc(cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),cwaveprj_idat,& 1050 & wi(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor,iblock),gs_hamk,mpi_enreg) 1051 mkl(1,fockcommon%ieigen,fockcommon%ieigen)=fockcommon%eigen_ikpt(fockcommon%ieigen) 1052 call timab(1563,2,tsec) ; call timab(1562,-1,tsec) 1053 fockcommon%e_fock0=fockcommon%e_fock0+half*weight(iblocksize)*fockcommon%eigen_ikpt(fockcommon%ieigen) 1054 if (fockcommon%optfor) then 1055 fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen) 1056 end if 1057 end do 1058 1059 1060 end do ! End of loop on block of bands 1061 1062 ! Calculate Mkl for the current k-point 1063 ABI_MALLOC(cwavefk,(2,npw_k*my_nspinor)) 1064 do iblock=1,nblockbd 1065 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 1066 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 1067 do iblocksize=1,blocksize 1068 kk=(iblock-1)*blocksize+iblocksize 1069 cwavefk(:,:)=cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor) 1070 do jblock=1,iblock 1071 do jblocksize=1,blocksize 1072 ll=(jblock-1)*blocksize+jblocksize 1073 if (ll<kk) then 1074 call dotprod_g(mkl(1,kk,ll),mkl(2,kk,ll),gs_hamk%istwf_k,npw_k,2,wi(:,1+(jblocksize-1)*npw_k*my_nspinor:& 1075 & jblocksize*npw_k*my_nspinor,jblock),cwavefk,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1076 end if 1077 end do 1078 end do 1079 end do 1080 end do ! End of loop on block of bands 1081 1082 ABI_FREE(cwavefk) 1083 mkl=-mkl 1084 1085 ! Cholesky factorisation of -mkl=Lx(trans(L)*. On output mkl=L 1086 call zpotrf("L",nband_k,mkl,nband_k,info) 1087 1088 ! calculate trans(L-1) 1089 ABI_MALLOC(bb,(2,nband_k,nband_k)) 1090 bb=zero 1091 do kk=1,nband_k 1092 bb(1,kk,kk)=one 1093 end do 1094 call ztrtrs("L","T","N",nband_k,nband_k,mkl,nband_k,bb,nband_k,info) 1095 fock%fockACE(ikpt,isppol)%xi=zero 1096 1097 ! Calculate ksi 1098 do kk=1,nband_k 1099 do jblock=1,nblockbd 1100 do jblocksize=1,blocksize 1101 ll=(jblock-1)*blocksize+jblocksize 1102 fock%fockACE(ikpt,isppol)%xi(1,:,kk)=fock%fockACE(ikpt,isppol)%xi(1,:,kk)+bb(1,ll,kk)*wi(1,1+(jblocksize-1)*& 1103 & npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)-& 1104 & bb(2,ll,kk)*wi(2,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock) 1105 fock%fockACE(ikpt,isppol)%xi(2,:,kk)=fock%fockACE(ikpt,isppol)%xi(2,:,kk)+bb(1,ll,kk)*wi(2,1+(jblocksize-1)*& 1106 npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)+& 1107 & bb(2,ll,kk)*wi(1,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock) 1108 end do 1109 end do 1110 end do 1111 1112 ! DEBUG 1113 ! fock%fockACE(ikpt,isppol)%xi=zero 1114 ! ENDDEBUG 1115 1116 ABI_FREE(wi) 1117 ABI_FREE(mkl) 1118 1119 ! Restore the bandfft tabs 1120 if (mpi_enreg%paral_kgb==1) then 1121 call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav) 1122 end if 1123 1124 ! Increment indices 1125 bdtot_index=bdtot_index+nband_k 1126 if (mkmem/=0) then 1127 ibg=ibg+my_nspinor*nband_cprj_k 1128 icg=icg+npw_k*my_nspinor*nband_k 1129 ikg=ikg+npw_k 1130 end if 1131 1132 if (psps%usepaw==1) then 1133 call pawcprj_free(cwaveprj) 1134 end if 1135 ABI_FREE(cwaveprj) 1136 ABI_FREE(cwavef) 1137 ABI_FREE(bb) 1138 ABI_FREE(occblock) 1139 ABI_FREE(weight) 1140 ABI_FREE(ffnl) 1141 ABI_FREE(kg_k) 1142 ABI_FREE(kpg_k) 1143 ABI_FREE(ylm_k) 1144 ABI_FREE(ylmgr_k) 1145 ABI_FREE(ph3d) 1146 end do ! End k point loop 1147 end do ! End loop over spins 1148 1149 call timab(1562,2,tsec) 1150 call timab(1565,1,tsec) 1151 1152 !Parallel case: accumulate (n,k) contributions 1153 if (xmpi_paral==1) then 1154 call xmpi_sum(fockcommon%e_fock0,spaceComm,ierr) 1155 ! Forces 1156 if (optfor==1) then 1157 if (psps%usepaw==1) then 1158 call xmpi_sum(fockcommon%forces,spaceComm,ierr) 1159 end if 1160 end if 1161 end if 1162 1163 !need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted) 1164 if (psps%usepaw==1) then 1165 call pawcprj_reorder(cprj,gs_hamk%atindx1) 1166 end if 1167 !Deallocate temporary space 1168 call gs_hamk%free() 1169 1170 call timab(1565,2,tsec) 1171 call timab(1560,2,tsec) 1172 1173 end subroutine fock2ACE
ABINIT/fock_ACE_getghc [ Functions ]
NAME
fock_ACE_getghc
FUNCTION
Compute the matrix elements <G|Vx|psi> of the Fock operator in the ACE context.
INPUTS
cwavef(2,npw*nspinor*ndat)= planewave coefficients of wavefunctions on which Fock operator is applied. gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied mpi_enreg= information about MPI parallelization
SIDE EFFECTS
ghc(2,npw*ndat)= matrix elements <G|H|C> or <G|H-lambda.S|C> (if sij_opt>=0 or =-1 in getghc) contains the fock exchange term for cwavef at the end.
NOTES
The current version assumes that : * nspinor = 1 * no "my_nspinor" * no restriction to the value of istwfk_bz (but must be tested in all case) * all the data for the occupied states (cgocc_bz) are the same as those for the current states (cg)
SOURCE
1201 subroutine fock_ACE_getghc(cwavef,ghc,gs_ham,mpi_enreg) 1202 1203 !Arguments ------------------------------------ 1204 ! Scalars 1205 type(MPI_type),intent(in) :: mpi_enreg 1206 type(gs_hamiltonian_type),target,intent(inout) :: gs_ham 1207 ! Arrays 1208 real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k) 1209 real(dp),intent(inout) :: ghc(:,:) 1210 1211 !Local variables------------------------------- 1212 ! Scalars 1213 integer :: iband,ikpt,ipw,my_nspinor,nband_k,npw 1214 real(dp) :: doti,dotr,eigen 1215 type(fock_common_type),pointer :: fockcommon 1216 ! Arrays 1217 real(dp) :: tsec(2) 1218 real(dp), allocatable :: ghc1(:,:),xi(:,:) 1219 1220 ! ************************************************************************* 1221 1222 call timab(1580,1,tsec) 1223 1224 ABI_CHECK(associated(gs_ham%fockcommon),"fock must be associated!") 1225 fockcommon => gs_ham%fockcommon 1226 1227 ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!") 1228 ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!") 1229 1230 ikpt=fockcommon%ikpt 1231 npw=gs_ham%npw_k 1232 nband_k=fockcommon%nband(ikpt) 1233 my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor) 1234 !*Initialization of the array ghc1 1235 !*ghc1 will contain the exact exchange contribution to the Hamiltonian 1236 ABI_MALLOC(ghc1,(2,npw*my_nspinor)) 1237 ghc1=zero 1238 ABI_MALLOC(xi,(2,npw*my_nspinor)) 1239 1240 do iband=1, nband_k 1241 xi(1,:)=gs_ham%fockACE_k%xi(1,:,iband) 1242 xi(2,:)=gs_ham%fockACE_k%xi(2,:,iband) 1243 1244 call dotprod_g(dotr,doti,gs_ham%istwf_k,npw*my_nspinor,2,xi,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1245 1246 ghc1(1,:)=ghc1(1,:)-(dotr*gs_ham%fockACE_k%xi(1,:,iband)-doti*gs_ham%fockACE_k%xi(2,:,iband)) 1247 ghc1(2,:)=ghc1(2,:)-(dotr*gs_ham%fockACE_k%xi(2,:,iband)+doti*gs_ham%fockACE_k%xi(1,:,iband)) 1248 end do 1249 ABI_FREE(xi) 1250 1251 !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc 1252 ! ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_spkpt + ghc1(:,:) 1253 ghc(:,:)=ghc(:,:) + ghc1(:,:) 1254 1255 ! call xmpi_sum(ghc,mpi_enreg%comm_kpt,ier) 1256 1257 ! ============================================ 1258 ! === Calculate the contribution to energy === 1259 ! ============================================ 1260 !* Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction 1261 !* at each convergence step. (consistent definition with the definition of hartree energy) 1262 if (fockcommon%ieigen/=0) then 1263 eigen=zero 1264 !* Dot product of cwavef and ghc 1265 !* inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering 1266 if(gs_ham%istwf_k==2) then 1267 eigen=half*cwavef(1,1)*ghc1(1,1) 1268 else 1269 eigen=cwavef(1,1)*ghc1(1,1)+cwavef(2,1)*ghc1(2,1) 1270 end if 1271 do ipw=2,npw 1272 eigen=eigen+cwavef(1,ipw)*ghc1(1,ipw)+cwavef(2,ipw)*ghc1(2,ipw) 1273 end do 1274 if(gs_ham%istwf_k>=2) eigen=two*eigen 1275 ! call xmpi_sum(eigen,mpi_enreg%comm_kpt,ier) 1276 fockcommon%eigen_ikpt(fockcommon%ieigen)= eigen 1277 fockcommon%ieigen = 0 1278 end if 1279 1280 ! =============================== 1281 ! === Deallocate local arrays === 1282 ! =============================== 1283 1284 ABI_FREE(ghc1) 1285 call timab(1580,2,tsec) 1286 1287 end subroutine fock_ACE_getghc
ABINIT/fock_getghc [ Functions ]
NAME
fock_getghc
FUNCTION
Compute the matrix elements <G|Vx|psi> of the Fock operator.
INPUTS
cwavef(2,npw*nspinor*ndat)= planewave coefficients of wavefunctions on which Fock operator is applied. cwaveprj <type(pawcprj_type> = <cwavevf|proj> gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied mpi_enreg= information about MPI parallelization
SIDE EFFECTS
ghc(2,npw*ndat)= matrix elements <G|H|C> or <G|H-lambda.S|C> (if sij_opt>=0 or =-1 in getghc) contains the fock exchange term for cwavef at the end.
NOTES
The current version assumes that: * nspinor = 1 * no "my_nspinor" * no restriction to the value of istwfk_bz (but must be tested in all case) * all the data for the occupied states (cgocc_bz) are the same as those for the current states (cg)
SOURCE
90 subroutine fock_getghc(cwavef,cwaveprj,ghc,gs_ham,mpi_enreg) 91 92 !Arguments ------------------------------------ 93 ! Scalars 94 type(MPI_type),intent(in) :: mpi_enreg 95 type(gs_hamiltonian_type),target,intent(inout) :: gs_ham 96 ! Arrays 97 type(pawcprj_type),intent(inout) :: cwaveprj(:,:) 98 real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k) 99 real(dp),intent(inout) :: ghc(:,:) 100 101 !Local variables------------------------------- 102 ! Scalars 103 integer,parameter :: tim_fourwf_fock_getghc=10,tim_fourdp_fock_getghc=10,ndat1=1 104 integer :: bdtot_jindex,choice,cplex_fock,cplex_dij,cpopt,i1,i2,i3,ia,iatom 105 integer :: iband_cprj,ider,idir,idir1,ier,ii,ind,ipw,ifft,itypat,izero,jband,jbg,jcg,jkg 106 integer :: jkpt,my_jsppol,jstwfk,lmn2_size,mgfftf,mpw,n1,n2,n3,n4,n5,n6 107 integer :: n1f,n2f,n3f,n4f,n5f,n6f,natom,nband_k,ndij,nfft,nfftf,nfftotf,nhat12_grdim,nnlout 108 integer :: npw,npwj,nspden_fock,nspinor,paw_opt,signs,tim_nonlop 109 integer, save :: ncount=0 110 logical :: need_ghc,qeq0 111 real(dp),parameter :: weight1=one 112 real(dp) :: doti,eigen,imcwf,imcwocc,imvloc,invucvol,recwf,recwocc,revloc,occ,wtk 113 type(fock_common_type),pointer :: fockcommon 114 type(fock_BZ_type),pointer :: fockbz 115 ! Arrays 116 integer :: ngfft(18),ngfftf(18) 117 integer,pointer :: gboundf(:,:),kg_occ(:,:),gbound_kp(:,:) 118 real(dp) :: enlout_dum(1),dotr(6),fockstr(6),for1(3),qphon(3),qvec_j(3),tsec(2),gsc_dum(2,0),rhodum(2,1) 119 real(dp) :: rhodum0(0,1,1),str(3,3) 120 real(dp), allocatable :: dummytab(:,:),dijhat(:,:,:,:),dijhat_tmp(:,:),ffnl_kp_dum(:,:,:,:) 121 real(dp), allocatable :: gvnlxc(:,:),ghc1(:,:),ghc2(:,:),grnhat12(:,:,:,:),grnhat_12(:,:,:,:,:),forikpt(:,:) 122 real(dp), allocatable :: rho12(:,:,:),rhog_munu(:,:),rhor_munu(:,:),vlocpsi_r(:) 123 real(dp), allocatable :: vfock(:),psilocal(:,:,:),vectin_dum(:,:),vqg(:),forout(:,:),strout(:,:) 124 real(dp), allocatable,target ::cwavef_r(:,:,:,:) 125 real(dp), ABI_CONTIGUOUS pointer :: cwaveocc_r(:,:,:,:) 126 type(pawcprj_type),pointer :: cwaveocc_prj(:,:) 127 128 real(dp) :: rprimd(3,3),for12(3) 129 130 ! ************************************************************************* 131 !return 132 133 ncount=ncount+1 134 135 call timab(1504,1,tsec) ; call timab(1505,-1,tsec) ; call timab(1515,-1,tsec) ; call timab(1541,-1,tsec) 136 137 ABI_CHECK(associated(gs_ham%fockcommon),"fock_common must be associated!") 138 fockcommon => gs_ham%fockcommon 139 ABI_CHECK(associated(gs_ham%fockbz),"fock_bz must be associated!") 140 fockbz => gs_ham%fockbz 141 142 ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!") 143 ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!") 144 if (fockcommon%usepaw==1) then 145 ABI_CHECK((size(cwaveprj,1)==gs_ham%natom.and.size(cwaveprj,2)==gs_ham%nspinor),"error on cwaveprj dims") 146 end if 147 need_ghc=(size(ghc,2)>0) 148 149 !Some constants 150 invucvol=1.d0/sqrt(gs_ham%ucvol) 151 call matr3inv(gs_ham%gprimd,rprimd) 152 cplex_fock=2;nspden_fock=1 153 natom=fockcommon%natom 154 nspinor=gs_ham%nspinor 155 mpw=maxval(fockbz%npwarr) 156 npw=gs_ham%npw_k 157 ider=0;izero=0 158 if (fockcommon%usepaw==1) then 159 nfft =fockcommon%pawfgr%nfftc ; ngfft =fockcommon%pawfgr%ngfftc 160 nfftf=fockcommon%pawfgr%nfft ; ngfftf=fockcommon%pawfgr%ngfft 161 mgfftf=fockcommon%pawfgr%mgfft 162 else 163 nfft =gs_ham%nfft ; nfftf =nfft 164 ngfft=gs_ham%ngfft ; ngfftf=ngfft 165 mgfftf=gs_ham%mgfft 166 end if 167 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 168 n4=ngfft(4);n5=ngfft(5);n6=ngfft(6) 169 n1f=ngfftf(1);n2f=ngfftf(2);n3f=ngfftf(3) 170 n4f=ngfftf(4);n5f=ngfftf(5);n6f=ngfftf(6) 171 172 ! =========================== 173 ! === Initialize arrays === 174 ! =========================== 175 ! transient optfor and optstress 176 ! fockcommon%optfor=.false. 177 ! fockcommon%optstr=.false. 178 !*Initialization of local pointers 179 !*Initialization of the array cwavef_r 180 !*cwavef_r = current wavefunction in r-space 181 ABI_MALLOC(cwavef_r,(2,n4f,n5f,n6f)) 182 !*rhormunu = overlap matrix between cwavef and (jkpt,mu) in R-space 183 ABI_MALLOC(rhor_munu,(cplex_fock,nfftf)) 184 !*rhogmunu = overlap matrix between cwavef and (jkpt,mu) in G-space 185 ABI_MALLOC(rhog_munu,(2,nfftf)) 186 !*dummytab = variables for fourwf 187 ABI_MALLOC(dummytab,(2,nfft)) 188 !*vfock = Fock potential 189 ABI_MALLOC(vfock,(cplex_fock*nfftf)) 190 !*vqg = 4pi/(G+q)**2 191 ABI_MALLOC(vqg,(nfftf)) 192 193 !*Initialization of the array ghc1 194 !*ghc1 will contain the exact exchange contribution to the Hamiltonian 195 ABI_MALLOC(ghc1,(2,npw)) 196 ABI_MALLOC(ghc2,(2,npw)) 197 ghc1=zero 198 ghc2=zero 199 !*Initialization of the array vlocpsi_r 200 !*vlocpsi_r = partial local Fock operator applied to cwavef in r-space and summed over all occupied (jkpt,mu) 201 ABI_MALLOC(vlocpsi_r,(cplex_fock*nfftf)) 202 vlocpsi_r=zero 203 204 !*Additional arrays in case of paw 205 if (fockcommon%usepaw==1) then 206 nhat12_grdim=0 207 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 208 ider=3 209 ABI_MALLOC(strout,(2,npw*nspinor)) 210 end if 211 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 212 ider=3 213 ABI_MALLOC(forout,(2,npw*nspinor)) 214 ABI_MALLOC(forikpt,(3,natom)) 215 forikpt=zero 216 end if 217 ABI_MALLOC(grnhat_12,(2,nfftf,nspinor**2,3,natom*(ider/3))) 218 ABI_MALLOC(gvnlxc,(2,npw*nspinor)) 219 ABI_MALLOC(grnhat12,(2,nfftf,nspinor**2,3*nhat12_grdim)) 220 end if 221 222 if (fockcommon%usepaw==1.or.fockcommon%optstr) then 223 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 224 call sphereboundary(gboundf,gs_ham%istwf_k,gs_ham%kg_k,mgfftf,npw) 225 else 226 gboundf=>gs_ham%gbound_k 227 end if 228 ! ========================================== 229 ! === Get cwavef in real space using FFT === 230 ! ========================================== 231 cwavef_r=zero 232 call timab(1515,2,tsec) ; call timab(1541,-2,tsec) ; call timab(1512,-1,tsec) 233 call fourwf(0,rhodum0,cwavef,rhodum,cwavef_r,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,& 234 & mgfftf,mpi_enreg,ndat1,ngfftf,npw,1,n4f,n5f,n6f,0,tim_fourwf_fock_getghc,weight1,weight1,& 235 & gpu_option=gs_ham%gpu_option) 236 call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1541,-1,tsec) 237 cwavef_r=cwavef_r*invucvol 238 239 ! ===================================================== 240 ! === Select the states in cgocc_bz with the same spin === 241 ! ===================================================== 242 !* Initialization of the indices/shifts, according to the value of isppol 243 !* bdtot_jindex = shift to be applied on the location of data in the array occ_bz ? 244 bdtot_jindex=0 245 !* jbg = shift to be applied on the location of data in the array cprj/occ 246 jbg=0;jcg=0 247 my_jsppol=fockcommon%isppol 248 if((fockcommon%isppol==2).and.(mpi_enreg%nproc_spkpt/=1)) my_jsppol=1 249 250 !=================================== 251 !=== Loop on the k-points in IBZ === 252 !=================================== 253 jkg=0 254 255 if (associated(gs_ham%ph3d_kp)) then 256 nullify (gs_ham%ph3d_kp) 257 end if 258 259 call timab(1505,2,tsec) ; call timab(1506,-1,tsec) ; call timab(1541,-2,tsec) 260 261 do jkpt=1,fockbz%mkpt 262 263 call timab(1521,1,tsec) 264 265 !* nband_k = number of bands at point k_j 266 nband_k=fockbz%nbandocc_bz(jkpt,my_jsppol) 267 !* wtk = weight in BZ of this k point 268 wtk=fockbz%wtk_bz(jkpt) !*sqrt(gs_ham%ucvol) 269 !* jstwfk= how is stored the wavefunction 270 jstwfk=fockbz%istwfk_bz(jkpt) 271 !* npwj= number of plane wave in basis for the wavefunction 272 npwj=fockbz%npwarr(jkpt) 273 !* Basis sphere of G vectors 274 if (allocated(fockbz%cgocc)) then 275 gbound_kp => fockbz%gbound_bz(:,:,jkpt) 276 kg_occ => fockbz%kg_bz(:,1+jkg:npwj+jkg) 277 end if 278 !* Load k^prime hamiltonian in the gs_ham datastructure 279 ! Note: ffnl_kp / ph3d_kp / gbound_kp are not used 280 281 if (associated(gs_ham%ph3d_kp)) then 282 ABI_MALLOC(gs_ham%ph3d_kp,(2,npwj,gs_ham%matblk)) 283 end if 284 285 call gs_ham%load_kprime(kpt_kp=fockbz%kptns_bz(:,jkpt),& 286 & istwf_kp=jstwfk,npw_kp=npwj,kg_kp=fockbz%kg_bz(:,1+jkg:npwj+jkg)) 287 !* Some temporary allocations needed for PAW 288 if (fockcommon%usepaw==1) then 289 ABI_MALLOC(vectin_dum,(2,npwj*nspinor)) 290 vectin_dum=zero 291 ABI_MALLOC(ffnl_kp_dum,(npwj,0,gs_ham%lmnmax,gs_ham%ntypat)) 292 call gs_ham%load_kprime(ffnl_kp=ffnl_kp_dum) 293 end if 294 295 ! ====================================== 296 ! === Calculate the vector q=k_i-k_j === 297 ! ====================================== 298 !* Evaluation of kpoint_j, the considered k-point in reduced coordinates 299 ! kpoint_j(:)=fockbz%kptns_bz(:,jkpt) 300 !* the vector qvec is expressed in reduced coordinates. 301 ! qvec(:)=kpoint_i(:)-kpoint_j(:) 302 qvec_j(:)=gs_ham%kpt_k(:)-fockbz%kptns_bz(:,jkpt) 303 qeq0=(qvec_j(1)**2+qvec_j(2)**2+qvec_j(3)**2<1.d-15) 304 call bare_vqg(qvec_j,fockcommon%gsqcut,gs_ham%gmet,fockcommon%usepaw,fockcommon%hyb_mixing,& 305 & fockcommon%hyb_mixing_sr,fockcommon%hyb_range_fock,nfftf,fockbz%nkpt_bz,ngfftf,gs_ham%ucvol,vqg) 306 307 call timab(1521,2,tsec) 308 309 ! ================================================= 310 ! === Loop on the band indices jband of cgocc_k === 311 ! ================================================= 312 do jband=1,nband_k 313 314 !* occ = occupancy of jband at this k point 315 occ=fockbz%occ_bz(jband+bdtot_jindex,my_jsppol) 316 if(occ<tol8) cycle 317 318 ! This timing is placed after the cycle ... 319 call timab(1522,1,tsec) ; call timab(1542,-1,tsec) 320 321 ! ============================================== 322 ! === Get cwaveocc_r in real space using FFT === 323 ! ============================================== 324 if (allocated(fockbz%cwaveocc_bz)) then 325 cwaveocc_r => fockbz%cwaveocc_bz(:,:,:,:,jband+jbg,my_jsppol) 326 else 327 ABI_MALLOC(cwaveocc_r,(2,n4f,n5f,n6f)) 328 cwaveocc_r=zero 329 call timab(1515,2,tsec) ; call timab(1512,-1,tsec) ; call timab(1542,-2,tsec) 330 call fourwf(1,rhodum0,fockbz%cgocc(:,1+jcg+npwj*(jband-1):jcg+jband*npwj,my_jsppol),rhodum,cwaveocc_r, & 331 & gbound_kp,gbound_kp,jstwfk,kg_occ,kg_occ,mgfftf,mpi_enreg,ndat1,ngfftf,& 332 & npwj,1,n4f,n5f,n6f,0,tim_fourwf_fock_getghc,weight1,weight1,gpu_option=gs_ham%gpu_option) 333 call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1542,-1,tsec) 334 cwaveocc_r=cwaveocc_r*invucvol 335 end if 336 337 ! ================================================ 338 ! === Get the overlap density matrix rhor_munu === 339 ! ================================================ 340 !* Calculate the overlap density matrix in real space = conj(cwaveocc_r)*cwavef_r 341 !* rhor_munu will contain the overlap density matrix. 342 ! vfock=-int{conj(cwaveocc_r)*cwavef_r*dr'/|r-r'|} 343 344 call timab(1522,2,tsec) ; call timab(1542,-2,tsec) ; call timab(1523,-1,tsec) 345 346 ind=0 347 do i3=1,n3f 348 do i2=1,n2f 349 do i1=1,n1f 350 ind=ind+1 351 recwf =cwavef_r(1,i1,i2,i3) ; imcwf =cwavef_r(2,i1,i2,i3) 352 recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3) 353 rhor_munu(1,ind)= recwocc*recwf+imcwocc*imcwf 354 rhor_munu(2,ind)= recwocc*imcwf-imcwocc*recwf 355 end do ! i1 356 end do ! i2 357 end do ! i3 358 359 call timab(1523,2,tsec) 360 361 ! ======================================================= 362 ! === Add compensation charge density in the PAW case === 363 ! ======================================================= 364 365 call timab(1524,-1,tsec) ; call timab(1544,-1,tsec) 366 367 if (fockcommon%usepaw==1) then 368 iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+jband 369 370 ABI_MALLOC(rho12,(2,nfftf,nspinor**2)) 371 372 cwaveocc_prj=>fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1) 373 374 call pawmknhat_psipsi(cwaveprj,cwaveocc_prj,ider,izero,natom,natom,nfftf,ngfftf,& 375 & nhat12_grdim,nspinor,fockcommon%ntypat,fockbz%pawang,fockcommon%pawfgrtab,grnhat12,rho12,& 376 & fockcommon%pawtab,gprimd=gs_ham%gprimd,grnhat_12=grnhat_12,qphon=qvec_j,xred=gs_ham%xred,atindx=gs_ham%atindx) 377 378 rhor_munu(1,:)=rhor_munu(1,:)+rho12(1,:,nspinor) 379 rhor_munu(2,:)=rhor_munu(2,:)-rho12(2,:,nspinor) 380 end if 381 call timab(1515,2,tsec) ; call timab(1513,-1,tsec) ; call timab(1544,-2,tsec) 382 ! Perform an FFT using fourwf to get rhog_munu = FFT^-1(rhor_munu) 383 call fourdp(cplex_fock,rhog_munu,rhor_munu,-1,mpi_enreg,nfftf,1,ngfftf,tim_fourdp_fock_getghc) 384 call timab(1513,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1544,-1,tsec) 385 386 if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 387 call strfock(gs_ham%gprimd,fockcommon%gsqcut,fockstr,fockcommon%hyb_mixing,fockcommon%hyb_mixing_sr,& 388 & fockcommon%hyb_range_fock,mpi_enreg,nfftf,ngfftf,fockbz%nkpt_bz,rhog_munu,gs_ham%ucvol,qvec_j) 389 fockcommon%stress_ikpt(:,fockcommon%ieigen)=fockcommon%stress_ikpt(:,fockcommon%ieigen)+fockstr(:)*occ*wtk 390 if (fockcommon%usepaw==0.and.(.not.need_ghc)) then 391 if (allocated(fockbz%cgocc)) then 392 ABI_FREE(cwaveocc_r) 393 end if 394 call timab(1524,2,tsec) ; call timab(1544,-2,tsec) 395 cycle 396 end if 397 end if 398 call timab(1524,2,tsec) ; call timab(1544,-2,tsec) 399 400 ! =================================================== 401 ! === Calculate the local potential vfockloc_munu === 402 ! =================================================== 403 !* Apply the Poisson solver to "rhog_munu" while taking into account the effect of the vector "qvec" 404 !* This is precisely what is done in the subroutine hartre, with option cplex=2. 405 !* vfock will contain the local Fock potential, the result of hartre routine. 406 !* vfock = FFT( rhog_munu/|g+qvec|^2 ) 407 call timab(1525,-1,tsec) ; call timab(1545,-1,tsec) 408 #if 0 409 410 call timab(1515,-2,tsec) ; call timab(1513,-1,tsec) 411 call hartre(cplex_fock,fockcommon%gsqcut,fockcommon%usepaw,mpi_enreg,nfftf,ngfftf,& 412 & mpi_enreg%paral_kgb,rhog_munu,rprimd,vfock,divgq0=fock%divgq0,qpt=qvec_j) 413 call timab(1513,2,tsec) ; call timab(1515,-1,tsec) 414 415 #else 416 do ifft=1,nfftf 417 rhog_munu(1,ifft) = rhog_munu(1,ifft) * vqg(ifft) 418 rhog_munu(2,ifft) = rhog_munu(2,ifft) * vqg(ifft) 419 end do 420 421 call timab(1515,2,tsec) ; call timab(1513,-1,tsec) ; call timab(1545,-2,tsec) 422 call fourdp(cplex_fock,rhog_munu,vfock,+1,mpi_enreg,nfftf,1,ngfftf,tim_fourdp_fock_getghc) 423 call timab(1513,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1545,-1,tsec) 424 #endif 425 call timab(1525,-2,tsec) ; call timab(1545,-2,tsec) 426 427 !=============================================================== 428 !======== Calculate Dij_Fock_hat contribution in case of PAW === 429 !=============================================================== 430 431 call timab(1526,-1,tsec) ; call timab(1546,-1,tsec) 432 433 if (fockcommon%usepaw==1) then 434 qphon=qvec_j;nfftotf=product(ngfftf(1:3)) 435 cplex_dij=1;ndij=nspden_fock 436 ABI_MALLOC(dijhat,(cplex_dij*gs_ham%dimekb1,natom,ndij,cplex_fock)) 437 dijhat=zero 438 do iatom=1,natom 439 itypat=gs_ham%typat(iatom) 440 lmn2_size=fockcommon%pawtab(itypat)%lmn2_size 441 ABI_MALLOC(dijhat_tmp,(cplex_fock*cplex_dij*lmn2_size,ndij)) 442 dijhat_tmp=zero 443 call pawdijhat(dijhat_tmp,cplex_dij,cplex_fock,gs_ham%gprimd,iatom,& 444 & natom,ndij,nfftf,nfftotf,nspden_fock,nspden_fock,fockbz%pawang,fockcommon%pawfgrtab(iatom),& 445 & fockcommon%pawtab(itypat),vfock,qphon,gs_ham%ucvol,gs_ham%xred) 446 do ii=1,cplex_fock 447 ind=(ii-1)*lmn2_size 448 dijhat(1:cplex_dij*lmn2_size,iatom,:,ii)=dijhat_tmp(ind+1:ind+cplex_dij*lmn2_size,:) 449 end do 450 ABI_FREE(dijhat_tmp) 451 end do 452 signs=2; cpopt=2;idir=0; paw_opt=1;nnlout=1;tim_nonlop=17 453 454 if(need_ghc) then 455 choice=1 456 call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec) 457 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 458 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,gvnlxc,enl=dijhat,& 459 & select_k=K_H_KPRIME) 460 call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec) 461 ghc2=ghc2-gvnlxc*occ*wtk 462 end if 463 464 ! Forces calculation 465 466 if (fockcommon%optfor.and.(fockcommon%ieigen/=0)) then 467 choice=2; dotr=zero;doti=zero;cpopt=4;tim_nonlop=17 468 do iatom=1,natom 469 do idir=1,3 470 call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec) 471 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 472 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,& 473 & forout,enl=dijhat,iatom_only=iatom,& 474 & select_k=K_H_KPRIME) 475 call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec) 476 call dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,forout,mpi_enreg%me_g0,mpi_enreg%comm_fft) 477 for1(idir)=zero 478 do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd 479 ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft) 480 for1(idir)=for1(idir)+vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-& 481 & vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom) 482 end do 483 end do 484 do idir=1,3 485 for12(idir)=rprimd(1,idir)*for1(1)+rprimd(2,idir)*for1(2)+rprimd(3,idir)*for1(3) 486 forikpt(idir,iatom)=forikpt(idir,iatom)-(for12(idir)*gs_ham%ucvol/nfftf+dotr(idir))*occ*wtk 487 end do 488 end do 489 end if 490 491 ! Stresses calculation 492 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 493 signs=2;choice=3;cpopt=4;tim_nonlop=17 494 495 ! first contribution 496 dotr=zero 497 do idir=1,6 498 call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec) 499 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 500 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,& 501 & strout,enl=dijhat,select_k=K_H_KPRIME) 502 call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec) 503 call dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,strout,mpi_enreg%me_g0,mpi_enreg%comm_fft) 504 fockcommon%stress_ikpt(idir,fockcommon%ieigen)=fockcommon%stress_ikpt(idir,fockcommon%ieigen)-& 505 & dotr(idir)*occ*wtk/gs_ham%ucvol 506 end do 507 ! second contribution 508 str=zero 509 do iatom=1,natom 510 do idir=1,3 511 do idir1=1,3 512 do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd 513 ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft) 514 str(idir,idir1)=str(idir,idir1)+(vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-& 515 & vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom))*fockcommon%pawfgrtab(iatom)%rfgd(idir1,ifft) 516 end do 517 end do 518 end do 519 end do 520 do idir=1,3 521 fockstr(idir)=str(idir,idir) 522 end do 523 fockstr(4)=(str(3,2)+str(2,3))*half 524 fockstr(5)=(str(3,1)+str(1,3))*half 525 fockstr(6)=(str(1,2)+str(2,1))*half 526 do idir=1,6 527 fockcommon%stress_ikpt(idir,fockcommon%ieigen)=fockcommon%stress_ikpt(idir,fockcommon%ieigen)+& 528 & fockstr(idir)/nfftf*occ*wtk 529 end do 530 531 ! third contribution 532 doti=zero 533 do ifft=1,nfftf 534 doti=doti+vfock(2*ifft-1)*rho12(1,ifft,nspinor)-vfock(2*ifft)*rho12(2,ifft,nspinor) 535 end do 536 fockcommon%stress_ikpt(1:3,fockcommon%ieigen)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen)-doti/nfftf*occ*wtk 537 ! doti=zero 538 ! do ifft=1,nfftf 539 ! doti=doti+vfock(2*ifft-1)*rhor_munu(1,ifft)-vfock(2*ifft)*rhor_munu(2,ifft) 540 ! end do 541 ! fockcommon%stress_ikpt(1:3,fockcommon%ieigen)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen)+doti/nfftf*occ*wtk*half 542 end if ! end stresses 543 544 ABI_FREE(dijhat) 545 ABI_FREE(rho12) 546 end if !end PAW 547 call timab(1526,2,tsec) ; call timab(1546,-2,tsec) 548 549 ! ============================================================= 550 ! === Apply the local potential vfockloc_munu to cwaveocc_r === 551 ! ============================================================= 552 call timab(1527,-1,tsec) 553 ind=0 554 do i3=1,ngfftf(3) 555 do i2=1,ngfftf(2) 556 do i1=1,ngfftf(1) 557 ind=ind+1 558 ! ind=i1+ngfftf(1)*(i2-1+ngfftf(2)*(i3-1)) 559 revloc=vfock(2*ind-1) ; imvloc=vfock(2*ind) 560 recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3) 561 vlocpsi_r(2*ind-1)=vlocpsi_r(2*ind-1)-(revloc*recwocc-imvloc*imcwocc)*occ*wtk 562 vlocpsi_r(2*ind )=vlocpsi_r(2*ind )-(revloc*imcwocc+imvloc*recwocc)*occ*wtk 563 end do 564 end do 565 end do 566 if (allocated(fockbz%cgocc)) then 567 ABI_FREE(cwaveocc_r) 568 end if 569 call timab(1527,2,tsec) 570 571 end do ! jband 572 573 ! ======================================================== 574 ! === End of loop : update of shifts and deallocations === 575 ! ==============================:========================= 576 !* Update of the shifts to be applied (reminder : mkmem is not 0, nspinor=1) 577 call timab(1528,1,tsec) 578 jcg=jcg+npwj*nband_k 579 jbg=jbg+nband_k 580 bdtot_jindex=bdtot_jindex+nband_k 581 jkg=jkg+npwj 582 if (fockcommon%usepaw==1) then 583 ABI_FREE(vectin_dum) 584 ABI_FREE(ffnl_kp_dum) 585 end if 586 if (associated(gs_ham%ph3d_kp)) then 587 ABI_FREE(gs_ham%ph3d_kp) 588 end if 589 call timab(1528,2,tsec) 590 591 end do ! jkpt 592 593 ! ======================================================== 594 ! === After loop === 595 ! ======================================================== 596 597 call timab(1506,2,tsec) ; call timab(1507,1,tsec) ; call timab(1547,-1,tsec) 598 599 if (fockcommon%usepaw==1) then 600 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 601 call timab(1547,2,tsec) ; call timab(1548,-1,tsec) 602 call xmpi_sum(forikpt,mpi_enreg%comm_hf,ier) 603 call timab(1548,2,tsec) ; call timab(1547,-1,tsec) 604 do iatom=1,natom !Loop over atom 605 ia=gs_ham%atindx(iatom) 606 fockcommon%forces_ikpt(:,ia,fockcommon%ieigen)=forikpt(:,iatom) 607 end do 608 end if 609 end if 610 if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 611 call timab(1547,2,tsec) ; call timab(1548,-1,tsec) 612 call xmpi_sum(fockcommon%stress_ikpt,mpi_enreg%comm_hf,ier) 613 call timab(1548,2,tsec) ; call timab(1547,-1,tsec) 614 end if 615 616 if (.not.need_ghc) then 617 618 ! =============================== 619 ! === Deallocate local arrays === 620 ! =============================== 621 ABI_FREE(cwavef_r) 622 ABI_FREE(ghc1) 623 ABI_FREE(ghc2) 624 ABI_FREE(rhor_munu) 625 ABI_FREE(rhog_munu) 626 ABI_FREE(vlocpsi_r) 627 ABI_FREE(dummytab) 628 ABI_FREE(vfock) 629 ABI_FREE(vqg) 630 if (fockcommon%usepaw==1) then 631 ABI_FREE(gvnlxc) 632 ABI_FREE(grnhat12) 633 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 634 ABI_FREE(forikpt) 635 ABI_FREE(forout) 636 end if 637 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 638 ABI_FREE(strout) 639 end if 640 ABI_FREE(grnhat_12) 641 end if 642 if(fockcommon%usepaw==1.or.fockcommon%optstr) then 643 ABI_FREE(gboundf) 644 end if 645 !*Restore gs_ham datastructure 646 647 if (associated(gs_ham%ph3d_kp)) then 648 ABI_MALLOC(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk)) 649 end if 650 call gs_ham%load_kprime(kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,& 651 & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k) 652 653 ! if (fockcommon%ieigen/=0) fockcommon%ieigen=0 654 655 else 656 657 ! *Restore gs_ham datastructure 658 659 if (associated(gs_ham%ph3d_kp)) then 660 ABI_MALLOC(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk)) 661 end if 662 call gs_ham%load_kprime(kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,& 663 & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k) 664 665 ! * Perform an FFT using fourwf to get ghc1 = FFT^-1(vlocpsi_r) 666 ABI_MALLOC(psilocal,(cplex_fock*n4f,n5f,n6f)) 667 call fftpac(1,mpi_enreg,nspden_fock,cplex_fock*n1f,n2f,n3f,cplex_fock*n4f,n5f,n6f,ngfft,vlocpsi_r,psilocal,2) 668 669 call timab(1515,2,tsec) ; call timab(1512,-1,tsec) ; call timab(1547,-2,tsec) 670 call fourwf(0,rhodum0,rhodum,ghc1,psilocal,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,& 671 & mgfftf,mpi_enreg,ndat1,ngfftf,1,npw,n4f,n5f,n6f,3,tim_fourwf_fock_getghc,weight1,weight1,& 672 & gpu_option=gs_ham%gpu_option) 673 call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1547,-1,tsec) 674 ABI_FREE(psilocal) 675 676 ghc1=ghc1*sqrt(gs_ham%ucvol)+ghc2 677 678 ! * If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc 679 ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_hf + ghc1(:,:) 680 681 call timab(1547,2,tsec) ; call timab(1548,-1,tsec) 682 call xmpi_sum(ghc,mpi_enreg%comm_hf,ier) 683 call timab(1548,2,tsec) ; call timab(1547,-1,tsec) 684 685 ! =============================== 686 ! === Deallocate local PAW arrays === 687 ! =============================== 688 if (fockcommon%usepaw==1) then 689 ABI_FREE(gvnlxc) 690 ABI_FREE(grnhat12) 691 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 692 ABI_FREE(forikpt) 693 ABI_FREE(forout) 694 end if 695 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 696 ABI_FREE(strout) 697 end if 698 ABI_FREE(grnhat_12) 699 end if 700 if(fockcommon%usepaw==1.or.fockcommon%optstr) then 701 ABI_FREE(gboundf) 702 end if 703 ! ============================================ 704 ! === Calculate the contribution to energy === 705 ! ============================================ 706 ! * Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction 707 ! * at each convergence step. (consistent definition with the definition of hartree energy) 708 if (fockcommon%ieigen/=0) then 709 eigen=zero 710 ! * Dot product of cwavef and ghc 711 ! * inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering 712 if(gs_ham%istwf_k==2) then 713 eigen=half*cwavef(1,1)*ghc1(1,1) 714 else 715 eigen=cwavef(1,1)*ghc1(1,1)+cwavef(2,1)*ghc1(2,1) 716 end if 717 do ipw=2,npw 718 eigen=eigen+cwavef(1,ipw)*ghc1(1,ipw)+cwavef(2,ipw)*ghc1(2,ipw) 719 end do 720 if(gs_ham%istwf_k>=2) eigen=two*eigen 721 call timab(1547,2,tsec) ; call timab(1548,-1,tsec) 722 call xmpi_sum(eigen,mpi_enreg%comm_hf,ier) 723 call timab(1548,2,tsec) ; call timab(1547,-1,tsec) 724 fockcommon%eigen_ikpt(fockcommon%ieigen)= eigen 725 if(fockcommon%use_ACE==0) fockcommon%ieigen = 0 726 end if 727 728 ! =============================== 729 ! === Deallocate local arrays === 730 ! =============================== 731 ABI_FREE(cwavef_r) 732 ABI_FREE(ghc1) 733 ABI_FREE(ghc2) 734 ABI_FREE(rhor_munu) 735 ABI_FREE(rhog_munu) 736 ABI_FREE(vlocpsi_r) 737 ABI_FREE(dummytab) 738 ABI_FREE(vfock) 739 ABI_FREE(vqg) 740 741 endif 742 743 call timab(1504,2,tsec) ; call timab(1507,-2,tsec) ; call timab(1515,-2,tsec) ; call timab(1547,-2,tsec) 744 745 end subroutine fock_getghc
ABINIT/m_fock_getghc [ Modules ]
NAME
FUNCTION
COPYRIGHT
Copyright (C) 2013-2024 ABINIT group (CMartins, FJ, MT, XG) 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
14 #if defined HAVE_CONFIG_H 15 #include "config.h" 16 #endif 17 18 #include "abi_common.h" 19 20 module m_fock_getghc 21 22 use defs_basis 23 use m_abicore 24 use m_errors 25 use m_xmpi 26 use m_fock 27 use m_pawcprj 28 !use m_cgtools 29 30 use defs_abitypes, only : mpi_type 31 use defs_datatypes, only : pseudopotential_type 32 use m_time, only : timab, time_accu 33 use m_symtk, only : matr3inv 34 use m_cgtools, only : dotprod_g 35 use m_kg, only : mkkpg 36 use m_fftcore, only : sphereboundary 37 use m_fft, only : fftpac, fourwf, fourdp 38 use m_fstrings, only : sjoin, itoa 39 use m_hamiltonian, only : gs_hamiltonian_type, K_H_KPRIME, init_hamiltonian 40 use m_pawdij, only : pawdijhat 41 use m_paw_nhat, only : pawmknhat_psipsi 42 use m_spacepar, only : hartre 43 use m_nonlop, only : nonlop 44 use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, bandfft_kpt_savetabs,bandfft_kpt_restoretabs, & 45 prep_bandfft_tabs 46 use m_pawtab, only : pawtab_type 47 use m_paw_ij, only : paw_ij_type 48 use m_mkffnl, only : mkffnl 49 use m_mpinfo, only : proc_distrb_cycle 50 51 implicit none 52 53 private