TABLE OF CONTENTS


ABINIT/fock2ACE [ Functions ]

[ Top ] [ 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
  use_gpu_cuda= 0 or 1 to know if we use cuda for nonlop call
  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

PARENTS

      scfcv

CHILDREN

      bandfft_kpt_restoretabs,bandfft_kpt_savetabs,destroy_hamiltonian
      dotprod_g,fock_getghc,init_hamiltonian,load_k_hamiltonian
      load_spin_hamiltonian,mkffnl,mkkpg,pawcprj_alloc,pawcprj_free
      pawcprj_get,pawcprj_reorder,prep_bandfft_tabs,timab,xmpi_sum,zpotrf
      ztrtrs

SOURCE

 801 subroutine fock2ACE(cg,cprj,fock,istwfk,kg,kpt,mband,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,&
 802 &  mpw,my_natom,natom,nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,&
 803 &  ntypat,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,typat,usecprj,use_gpu_cuda,wtk,xred,ylm)
 804 
 805 
 806 !This section has been created automatically by the script Abilint (TD).
 807 !Do not modify the following lines by hand.
 808 #undef ABI_FUNC
 809 #define ABI_FUNC 'fock2ACE'
 810 !End of the abilint section
 811 
 812  implicit none
 813 
 814 !Arguments ------------------------------------
 815 !scalars
 816  integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nkpt
 817  integer,intent(in) :: nspden,nsppol,nspinor,ntypat,optfor
 818  integer,intent(in) :: usecprj,use_gpu_cuda
 819  type(MPI_type),intent(inout) :: mpi_enreg
 820  type(pseudopotential_type),intent(in) :: psps
 821 !arrays
 822  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
 823  integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt)
 824  integer,intent(in) :: typat(natom)
 825  real(dp),intent(in) :: cg(2,mcg)
 826  real(dp),intent(in) :: kpt(3,nkpt)
 827  real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
 828  real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom)
 829  real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
 830  type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
 831  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
 832  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
 833  type(fock_type),pointer, intent(inout) :: fock
 834 !Local variables-------------------------------
 835 !scalars
 836  integer :: bandpp,bdtot_index,dimffnl,iband,iband_cprj,iband_last,ibg,icg,ider
 837  integer :: idir,ierr,ikg,ikpt,ilm,ipw,isppol,istwf_k,kk,ll
 838  integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,ndat,nkpg
 839  integer :: npw_k,spaceComm
 840  integer :: use_ACE_old
 841  integer :: blocksize,iblock,jblock,iblocksize,jblocksize,nblockbd
 842 !integer, save :: counter=0
 843  type(gs_hamiltonian_type) :: gs_hamk
 844  logical :: compute_gbound
 845  character(len=500) :: msg
 846  type(fock_common_type),pointer :: fockcommon
 847 !arrays
 848  integer,allocatable :: kg_k(:,:)
 849  real(dp) :: kpoint(3),rmet(3,3),tsec(2)
 850  real(dp),allocatable :: bb(:,:,:),cwavef(:,:),cwavefk(:,:),ffnl_sav(:,:,:,:)
 851  real(dp),allocatable :: kpg_k(:,:),kpg_k_sav(:,:)
 852  real(dp),allocatable :: mkl(:,:,:),occblock(:),ph3d(:,:,:),ph3d_sav(:,:,:)
 853  real(dp),allocatable :: wi(:,:,:),weight(:),ylm_k(:,:),ylmgr_k(:,:,:)
 854  real(dp),allocatable,target :: ffnl(:,:,:,:)
 855  type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
 856  type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
 857  type(pawcprj_type),pointer :: cwaveprj_idat(:,:)
 858 
 859 !*************************************************************************
 860 
 861  call timab(920,1,tsec)
 862  call timab(921,1,tsec)
 863 
 864 !DEBUG
 865 !if(counter>0)return
 866 !counter=counter+1
 867 !ENDDEBUG
 868 
 869 !Init mpicomm and me
 870  if(mpi_enreg%paral_kgb==1)then
 871    spaceComm=mpi_enreg%comm_kpt
 872    me_distrb=mpi_enreg%me_kpt
 873  else
 874 !* In case of HF calculation
 875    if (mpi_enreg%paral_hf==1) then
 876      spaceComm=mpi_enreg%comm_kpt
 877      me_distrb=mpi_enreg%me_kpt
 878    else
 879      spaceComm=mpi_enreg%comm_cell
 880      me_distrb=mpi_enreg%me_cell
 881    end if
 882  end if
 883 
 884 !Some initializations
 885  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
 886  compute_gbound=.true.
 887  fockcommon => fock%fock_common
 888  use_ACE_old=fockcommon%use_ACE
 889  fockcommon%use_ACE=0
 890 
 891 !Initialize Hamiltonian (k- and spin-independent terms)
 892 
 893  call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,&
 894 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj,&
 895 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 896 & paw_ij=paw_ij,ph1d=ph1d,fock=fock,&
 897 & use_gpu_cuda=use_gpu_cuda)
 898  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
 899  fockcommon%use_ACE=use_ACE_old
 900  call timab(921,2,tsec)
 901 
 902 !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted)
 903  if (psps%usepaw==1) then
 904    call pawcprj_reorder(cprj,gs_hamk%atindx)
 905  end if
 906 
 907 !LOOP OVER SPINS
 908  bdtot_index=0;ibg=0;icg=0
 909 
 910  do isppol=1,nsppol
 911    fockcommon%isppol=isppol
 912 !  Continue to initialize the Hamiltonian (PAW DIJ coefficients)
 913    call load_spin_hamiltonian(gs_hamk,isppol,with_nonlocal=.true.)
 914 
 915 !  Loop over k points
 916    ikg=0
 917    do ikpt=1,nkpt
 918      fockcommon%ikpt=ikpt
 919      nband_k=nband(ikpt+(isppol-1)*nkpt)
 920      npw_k=npwarr(ikpt)
 921      kpoint(:)=kpt(:,ikpt)
 922      istwf_k=istwfk(ikpt)
 923      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
 924        bdtot_index=bdtot_index+nband_k
 925        cycle
 926      end if
 927 
 928      call timab(922,1,tsec)
 929 
 930 !    Parallelism over FFT and/or bands: define sizes and tabs
 931      if (mpi_enreg%paral_kgb==1) then
 932        my_ikpt=mpi_enreg%my_kpttab(ikpt)
 933        nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
 934        bandpp=mpi_enreg%bandpp
 935        my_bandfft_kpt => bandfft_kpt(my_ikpt)
 936      else
 937        my_ikpt=ikpt
 938        bandpp=mpi_enreg%bandpp
 939        nblockbd=nband_k/bandpp
 940      end if
 941      blocksize=nband_k/nblockbd
 942      mband_cprj=mband/mpi_enreg%nproc_band
 943      nband_cprj_k=nband_k/mpi_enreg%nproc_band
 944 
 945      ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize))
 946      if (psps%usepaw==1) then
 947        ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,my_nspinor*bandpp))
 948        call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
 949      else
 950        ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0))
 951      end if
 952 
 953      ABI_ALLOCATE(kg_k,(3,mpw))
 954 !$OMP PARALLEL DO
 955      do ipw=1,npw_k
 956        kg_k(:,ipw)=kg(:,ipw+ikg)
 957      end do
 958 
 959      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
 960      ABI_ALLOCATE(ylmgr_k,(0,0,0))
 961      if (psps%useylm==1) then
 962 !$OMP PARALLEL DO COLLAPSE(2)
 963        do ilm=1,mpsang*mpsang
 964          do ipw=1,npw_k
 965            ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
 966          end do
 967        end do
 968      end if
 969 
 970 !    Compute (k+G) vectors
 971      nkpg=3*nloalg(3)
 972      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 973      if (nkpg>0) then
 974        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 975      end if
 976 
 977 
 978 !    Compute nonlocal form factors ffnl at all (k+G)
 979      ider=0;idir=0;dimffnl=1
 980      ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
 981      call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
 982 &     ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
 983 &     nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
 984 
 985 !    Load k-dependent part in the Hamiltonian datastructure
 986 !     - Compute 3D phase factors
 987 !     - Prepare various tabs in case of band-FFT parallelism
 988 !     - Load k-dependent quantities in the Hamiltonian
 989 
 990      ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
 991      call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
 992 &     kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.)
 993 
 994 !    Load band-FFT tabs (transposed k-dependent arrays)
 995      if (mpi_enreg%paral_kgb==1) then
 996        call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
 997        call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg)
 998        call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, &
 999 &       kg_k     =my_bandfft_kpt%kg_k_gather, &
1000 &       kpg_k    =my_bandfft_kpt%kpg_k_gather, &
1001        ffnl_k   =my_bandfft_kpt%ffnl_gather, &
1002        ph3d_k   =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound)
1003      end if
1004 
1005      call timab(922,2,tsec)
1006 
1007 !    The following is now wrong. In sequential, nblockbd=nband_k/bandpp
1008 !    blocksize= bandpp (JB 2016/04/16)
1009 !    Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
1010 !
1011      ABI_ALLOCATE(occblock,(blocksize))
1012      ABI_ALLOCATE(weight,(blocksize))
1013      occblock=zero;weight=zero
1014 
1015      if (fockcommon%optfor) then
1016        fockcommon%forces_ikpt=zero
1017      end if
1018 
1019      ABI_ALLOCATE(wi,(2,npw_k*my_nspinor*blocksize,nblockbd))
1020      wi=zero
1021      ABI_ALLOCATE(mkl,(2,nband_k,nband_k))
1022      mkl=zero
1023 ! Calculate all the Wi for the current k-point
1024 
1025      do iblock=1,nblockbd
1026 
1027        iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
1028        iband_cprj=(iblock-1)*bandpp+1
1029        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
1030 
1031 !      Select occupied bandsddk
1032        occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1033        call timab(923,1,tsec)
1034        weight(:)=wtk(ikpt)*occblock(:)
1035 
1036 !        Load contribution from n,k
1037        cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1038 &       cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1039        if (psps%usepaw==1) then
1040          call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,&
1041 &         mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,&
1042 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1043        end if
1044 
1045        call timab(926,2,tsec)
1046 
1047        if (mpi_enreg%paral_kgb==1) then
1048          msg='fock2ACE: Paral_kgb is not yet implemented for fock calculations'
1049          MSG_BUG(msg)
1050        end if
1051        ndat=mpi_enreg%bandpp
1052        if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj
1053 
1054        do iblocksize=1,blocksize
1055          fockcommon%ieigen=(iblock-1)*blocksize+iblocksize
1056          fockcommon%iband=(iblock-1)*blocksize+iblocksize
1057          if (gs_hamk%usepaw==1) then
1058            cwaveprj_idat => cwaveprj(:,(iblocksize-1)*my_nspinor+1:iblocksize*my_nspinor)
1059          end if
1060          call fock_getghc(cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),cwaveprj_idat,&
1061 &         wi(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor,iblock),gs_hamk,mpi_enreg)
1062          mkl(1,fockcommon%ieigen,fockcommon%ieigen)=fockcommon%eigen_ikpt(fockcommon%ieigen)
1063          if (fockcommon%optfor) then
1064            fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen)
1065          end if
1066        end do
1067 
1068 
1069      end do ! End of loop on block of bands
1070 
1071 ! Calculate Mkl for the current k-point
1072      ABI_ALLOCATE(cwavefk,(2,npw_k*my_nspinor))
1073      do iblock=1,nblockbd
1074        cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1075 &       cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1076        do iblocksize=1,blocksize
1077          kk=(iblock-1)*blocksize+iblocksize
1078          cwavefk(:,:)=cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor)
1079          do jblock=1,iblock
1080            do jblocksize=1,blocksize
1081              ll=(jblock-1)*blocksize+jblocksize
1082              if (ll<kk) then
1083                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:&
1084 &               jblocksize*npw_k*my_nspinor,jblock),cwavefk,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1085              end if
1086            end do
1087          end do
1088        end do
1089      end do ! End of loop on block of bands
1090 
1091      ABI_DEALLOCATE(cwavefk)
1092      mkl=-mkl
1093 
1094 ! Cholesky factorisation of -mkl=Lx(trans(L)*. On output mkl=L
1095      call zpotrf("L",nband_k,mkl,nband_k,ierr)
1096 
1097 ! calculate trans(L-1)
1098      ABI_ALLOCATE(bb,(2,nband_k,nband_k))
1099      bb=zero
1100      do kk=1,nband_k
1101        bb(1,kk,kk)=one
1102      end do
1103      call ztrtrs("L","T","N",nband_k,nband_k,mkl,nband_k,bb,nband_k,ierr)
1104      fock%fockACE(ikpt,isppol)%xi=zero
1105 
1106 ! Calculate ksi
1107      do kk=1,nband_k
1108        do jblock=1,nblockbd
1109          do jblocksize=1,blocksize
1110            ll=(jblock-1)*blocksize+jblocksize
1111            fock%fockACE(ikpt,isppol)%xi(1,:,kk)=fock%fockACE(ikpt,isppol)%xi(1,:,kk)+bb(1,ll,kk)*wi(1,1+(jblocksize-1)*&
1112 &           npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)-&
1113 &           bb(2,ll,kk)*wi(2,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)
1114            fock%fockACE(ikpt,isppol)%xi(2,:,kk)=fock%fockACE(ikpt,isppol)%xi(2,:,kk)+bb(1,ll,kk)*wi(2,1+(jblocksize-1)*&
1115            npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)+&
1116 &           bb(2,ll,kk)*wi(1,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)
1117          end do
1118        end do
1119      end do
1120 
1121 !    DEBUG
1122 !    fock%fockACE(ikpt,isppol)%xi=zero
1123 !    ENDDEBUG
1124 
1125      ABI_DEALLOCATE(wi)
1126      ABI_DEALLOCATE(mkl)
1127 
1128 !    Restore the bandfft tabs
1129      if (mpi_enreg%paral_kgb==1) then
1130        call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1131      end if
1132 
1133 !    Increment indices
1134      bdtot_index=bdtot_index+nband_k
1135      if (mkmem/=0) then
1136        ibg=ibg+my_nspinor*nband_cprj_k
1137        icg=icg+npw_k*my_nspinor*nband_k
1138        ikg=ikg+npw_k
1139      end if
1140 
1141      if (psps%usepaw==1) then
1142        call pawcprj_free(cwaveprj)
1143      end if
1144      ABI_DATATYPE_DEALLOCATE(cwaveprj)
1145      ABI_DEALLOCATE(cwavef)
1146      ABI_DEALLOCATE(bb)
1147      ABI_DEALLOCATE(occblock)
1148      ABI_DEALLOCATE(weight)
1149      ABI_DEALLOCATE(ffnl)
1150      ABI_DEALLOCATE(kg_k)
1151      ABI_DEALLOCATE(kpg_k)
1152      ABI_DEALLOCATE(ylm_k)
1153      ABI_DEALLOCATE(ylmgr_k)
1154      ABI_DEALLOCATE(ph3d)
1155    end do ! End k point loop
1156  end do ! End loop over spins
1157 
1158 !Parallel case: accumulate (n,k) contributions
1159  if (xmpi_paral==1) then
1160 !  Forces
1161    if (optfor==1) then
1162      call timab(65,2,tsec)
1163      if (psps%usepaw==1) then
1164        call xmpi_sum(fockcommon%forces,spaceComm,ierr)
1165      end if
1166    end if
1167  end if
1168 
1169  call timab(925,1,tsec)
1170 
1171 !need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted)
1172  if (psps%usepaw==1) then
1173    call pawcprj_reorder(cprj,gs_hamk%atindx1)
1174  end if
1175 !Deallocate temporary space
1176  call destroy_hamiltonian(gs_hamk)
1177 
1178  call timab(925,2,tsec)
1179  call timab(920,2,tsec)
1180 
1181 end subroutine fock2ACE

ABINIT/fock_ACE_getghc [ Functions ]

[ Top ] [ 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)

PARENTS

      getghc

CHILDREN

      dotprod_g

SOURCE

1215 subroutine fock_ACE_getghc(cwavef,ghc,gs_ham,mpi_enreg)
1216 
1217 
1218 !This section has been created automatically by the script Abilint (TD).
1219 !Do not modify the following lines by hand.
1220 #undef ABI_FUNC
1221 #define ABI_FUNC 'fock_ACE_getghc'
1222 !End of the abilint section
1223 
1224  implicit none
1225 
1226 !Arguments ------------------------------------
1227 ! Scalars
1228  type(MPI_type),intent(in) :: mpi_enreg
1229  type(gs_hamiltonian_type),target,intent(inout) :: gs_ham
1230 ! Arrays
1231  real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k)
1232  real(dp),intent(inout) :: ghc(:,:)
1233 
1234 !Local variables-------------------------------
1235 ! Scalars
1236  integer :: iband,ikpt,ipw,my_nspinor,nband_k,npw
1237  real(dp) :: doti,dotr,eigen
1238  type(fock_common_type),pointer :: fockcommon
1239 ! Arrays
1240  real(dp), allocatable :: ghc1(:,:),xi(:,:)
1241 
1242 
1243 ! *************************************************************************
1244 
1245  ABI_CHECK(associated(gs_ham%fockcommon),"fock must be associated!")
1246  fockcommon => gs_ham%fockcommon
1247 
1248  ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!")
1249  ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!")
1250 
1251  ikpt=fockcommon%ikpt
1252  npw=gs_ham%npw_k
1253  nband_k=fockcommon%nband(ikpt)
1254  my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor)
1255 !*Initialization of the array ghc1
1256 !*ghc1 will contain the exact exchange contribution to the Hamiltonian
1257  ABI_ALLOCATE(ghc1,(2,npw*my_nspinor))
1258  ghc1=zero
1259  ABI_ALLOCATE(xi,(2,npw*my_nspinor))
1260 
1261  do iband=1, nband_k
1262    xi(1,:)=gs_ham%fockACE_k%xi(1,:,iband)
1263    xi(2,:)=gs_ham%fockACE_k%xi(2,:,iband)
1264 
1265    call dotprod_g(dotr,doti,gs_ham%istwf_k,npw*my_nspinor,2,xi,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1266 
1267    ghc1(1,:)=ghc1(1,:)-(dotr*gs_ham%fockACE_k%xi(1,:,iband)-doti*gs_ham%fockACE_k%xi(2,:,iband))
1268    ghc1(2,:)=ghc1(2,:)-(dotr*gs_ham%fockACE_k%xi(2,:,iband)+doti*gs_ham%fockACE_k%xi(1,:,iband))
1269  end do
1270  ABI_DEALLOCATE(xi)
1271 
1272 !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc
1273 ! ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_kpt + ghc1(:,:)
1274  ghc(:,:)=ghc(:,:) + ghc1(:,:)
1275 
1276 ! call xmpi_sum(ghc,mpi_enreg%comm_kpt,ier)
1277 
1278 ! ============================================
1279 ! === Calculate the contribution to energy ===
1280 ! ============================================
1281 !* Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction
1282 !* at each convergence step. (consistent definition with the defintion of hartree energy)
1283  if (fockcommon%ieigen/=0) then
1284    eigen=zero
1285 !* Dot product of cwavef and ghc
1286 !* inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering
1287    if(gs_ham%istwf_k==2) then
1288      eigen=half*cwavef(1,1)*ghc1(1,1)
1289    else
1290      eigen=cwavef(1,1)*ghc1(1,1)+cwavef(2,1)*ghc1(2,1)
1291    end if
1292    do ipw=2,npw
1293      eigen=eigen+cwavef(1,ipw)*ghc1(1,ipw)+cwavef(2,ipw)*ghc1(2,ipw)
1294    end do
1295    if(gs_ham%istwf_k>=2) eigen=two*eigen
1296 !   call xmpi_sum(eigen,mpi_enreg%comm_kpt,ier)
1297    fockcommon%eigen_ikpt(fockcommon%ieigen)= eigen
1298    fockcommon%ieigen = 0
1299  end if
1300 
1301 ! ===============================
1302 ! === Deallocate local arrays ===
1303 ! ===============================
1304 
1305  ABI_DEALLOCATE(ghc1)
1306 
1307 end subroutine fock_ACE_getghc

ABINIT/fock_getghc [ Functions ]

[ Top ] [ 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)

PARENTS

      fock2ACE,forstrnps,getghc

CHILDREN

      bare_vqg,dotprod_g,fftpac,fourdp,fourwf,hartre,load_kprime_hamiltonian
      matr3inv,nonlop,pawdijhat,pawmknhat_psipsi,sphereboundary,strfock,timab
      xmpi_sum

SOURCE

104 subroutine fock_getghc(cwavef,cwaveprj,ghc,gs_ham,mpi_enreg)
105 
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'fock_getghc'
111 !End of the abilint section
112 
113  implicit none
114 
115 !Arguments ------------------------------------
116 ! Scalars
117  type(MPI_type),intent(in) :: mpi_enreg
118  type(gs_hamiltonian_type),target,intent(inout) :: gs_ham
119 ! Arrays
120  type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
121  real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k)
122  real(dp),intent(inout) :: ghc(:,:)
123 
124 !Local variables-------------------------------
125 ! Scalars
126  integer,parameter :: tim_fourwf0=0,tim_fourdp0=0,ndat1=1
127  integer :: bdtot_jindex,choice,cplex_fock,cplex_dij,cpopt,i1,i2,i3,ia,iatom
128  integer :: iband_cprj,ider,idir,idir1,ier,ii,ind,ipert,ipw,ifft,itypat,izero,jband,jbg,jcg,jkg
129  integer :: jkpt,my_jsppol,jstwfk,lmn2_size,mgfftf,mpw,n1,n2,n3,n4,n5,n6
130  integer :: n1f,n2f,n3f,n4f,n5f,n6f,natom,nband_k,ndij,nfft,nfftf,nfftotf,nhat12_grdim,nnlout
131  integer :: npw,npwj,nspden_fock,nspinor,paw_opt,signs,tim_nonlop
132  logical :: need_ghc,qeq0
133  real(dp),parameter :: weight1=one
134  real(dp) :: doti,eigen,imcwf,imcwocc,imvloc,invucvol,recwf,recwocc,revloc,occ,wtk
135  type(fock_common_type),pointer :: fockcommon
136  type(fock_BZ_type),pointer :: fockbz
137 ! Arrays
138  integer :: ngfft(18),ngfftf(18)
139  integer,pointer :: gboundf(:,:),kg_occ(:,:),gbound_kp(:,:)
140  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)
141  real(dp) :: rhodum0(0,1,1),str(3,3)
142  real(dp), allocatable :: dummytab(:,:),dijhat(:,:,:,:),dijhat_tmp(:,:),ffnl_kp_dum(:,:,:,:)
143  real(dp), allocatable :: gvnlc(:,:),ghc1(:,:),ghc2(:,:),grnhat12(:,:,:,:),grnhat_12(:,:,:,:,:),forikpt(:,:)
144  real(dp), allocatable :: rho12(:,:,:),rhog_munu(:,:),rhor_munu(:,:),vlocpsi_r(:)
145  real(dp), allocatable :: vfock(:),psilocal(:,:,:),vectin_dum(:,:),vqg(:),forout(:,:),strout(:,:)
146  real(dp), allocatable,target ::cwavef_r(:,:,:,:)
147  real(dp), ABI_CONTIGUOUS  pointer :: cwaveocc_r(:,:,:,:)
148  type(pawcprj_type),pointer :: cwaveocc_prj(:,:)
149 
150  real(dp) :: rprimd(3,3),for12(3)
151 
152 ! *************************************************************************
153 !return
154  call timab(1504,1,tsec)
155  call timab(1505,1,tsec)
156 
157  ABI_CHECK(associated(gs_ham%fockcommon),"fock_common must be associated!")
158  fockcommon => gs_ham%fockcommon
159  ABI_CHECK(associated(gs_ham%fockbz),"fock_bz must be associated!")
160  fockbz => gs_ham%fockbz
161 
162  ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!")
163  ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!")
164  if (fockcommon%usepaw==1) then
165    ABI_CHECK((size(cwaveprj,1)==gs_ham%natom.and.size(cwaveprj,2)==gs_ham%nspinor),"error on cwaveprj dims")
166  end if
167  need_ghc=(size(ghc,2)>0)
168 
169 !Some constants
170  invucvol=1.d0/sqrt(gs_ham%ucvol)
171  call matr3inv(gs_ham%gprimd,rprimd)
172  cplex_fock=2;nspden_fock=1
173  natom=fockcommon%natom
174  nspinor=gs_ham%nspinor
175  mpw=maxval(fockbz%npwarr)
176  npw=gs_ham%npw_k
177  ider=0;izero=0
178  if (fockcommon%usepaw==1) then
179    nfft =fockcommon%pawfgr%nfftc ; ngfft =fockcommon%pawfgr%ngfftc
180    nfftf=fockcommon%pawfgr%nfft  ; ngfftf=fockcommon%pawfgr%ngfft
181    mgfftf=fockcommon%pawfgr%mgfft
182  else
183    nfft =gs_ham%nfft  ; nfftf =nfft
184    ngfft=gs_ham%ngfft ; ngfftf=ngfft
185    mgfftf=gs_ham%mgfft
186  end if
187  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
188  n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
189  n1f=ngfftf(1);n2f=ngfftf(2);n3f=ngfftf(3)
190  n4f=ngfftf(4);n5f=ngfftf(5);n6f=ngfftf(6)
191 
192 ! ===========================
193 ! === Initialize arrays   ===
194 ! ===========================
195 ! transient optfor and optstress
196 ! fockcommon%optfor=.false.
197 ! fockcommon%optstr=.false.
198 !*Initialization of local pointers
199 !*Initialization of the array cwavef_r
200 !*cwavef_r = current wavefunction in r-space
201  ABI_ALLOCATE(cwavef_r,(2,n4f,n5f,n6f))
202 !*rhormunu = overlap matrix between cwavef and (jkpt,mu) in R-space
203  ABI_ALLOCATE(rhor_munu,(cplex_fock,nfftf))
204 !*rhogmunu = overlap matrix between cwavef and (jkpt,mu) in G-space
205  ABI_ALLOCATE(rhog_munu,(2,nfftf))
206 !*dummytab = variables for fourwf
207  ABI_ALLOCATE(dummytab,(2,nfft))
208 !*vfock = Fock potential
209  ABI_ALLOCATE(vfock,(cplex_fock*nfftf))
210 !*vqg = 4pi/(G+q)**2
211  ABI_ALLOCATE(vqg,(nfftf))
212 
213 !*Initialization of the array ghc1
214 !*ghc1 will contain the exact exchange contribution to the Hamiltonian
215  ABI_ALLOCATE(ghc1,(2,npw))
216  ABI_ALLOCATE(ghc2,(2,npw))
217  ghc1=zero
218  ghc2=zero
219 !*Initialization of the array vlocpsi_r
220 !*vlocpsi_r = partial local Fock operator applied to cwavef in r-space and summed over all occupied (jkpt,mu)
221  ABI_ALLOCATE(vlocpsi_r,(cplex_fock*nfftf))
222  vlocpsi_r=zero
223 
224 !*Additional arrays in case of paw
225  if (fockcommon%usepaw==1) then
226    nhat12_grdim=0
227    if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
228      ider=3
229      ABI_ALLOCATE(strout,(2,npw*nspinor))
230    end if
231    if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
232      ider=3
233      ABI_ALLOCATE(forout,(2,npw*nspinor))
234      ABI_ALLOCATE(forikpt,(3,natom))
235      forikpt=zero
236    end if
237    ABI_ALLOCATE(grnhat_12,(2,nfftf,nspinor**2,3,natom*(ider/3)))
238    ABI_ALLOCATE(gvnlc,(2,npw*nspinor))
239    ABI_ALLOCATE(grnhat12,(2,nfftf,nspinor**2,3*nhat12_grdim))
240  end if
241 
242  if (fockcommon%usepaw==1.or.fockcommon%optstr) then
243    ABI_ALLOCATE(gboundf,(2*mgfftf+8,2))
244    call sphereboundary(gboundf,gs_ham%istwf_k,gs_ham%kg_k,mgfftf,npw)
245  else
246    gboundf=>gs_ham%gbound_k
247  end if
248 
249 ! ==========================================
250 ! === Get cwavef in real space using FFT ===
251 ! ==========================================
252  cwavef_r=zero
253  call fourwf(0,rhodum0,cwavef,rhodum,cwavef_r,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,&
254 & mgfftf,mpi_enreg,ndat1,ngfftf,npw,1,n4f,n5f,n6f,0,mpi_enreg%paral_kgb,tim_fourwf0,weight1,weight1,&
255 & use_gpu_cuda=gs_ham%use_gpu_cuda)
256  cwavef_r=cwavef_r*invucvol
257 
258 ! =====================================================
259 ! === Select the states in cgocc_bz with the same spin ===
260 ! =====================================================
261 !* Initialization of the indices/shifts, according to the value of isppol
262 !* bdtot_jindex = shift to be applied on the location of data in the array occ_bz ?
263  bdtot_jindex=0
264 !* jbg = shift to be applied on the location of data in the array cprj/occ
265  jbg=0;jcg=0
266  my_jsppol=fockcommon%isppol
267  if((fockcommon%isppol==2).and.(mpi_enreg%nproc_kpt/=1)) my_jsppol=1
268 
269  call timab(1505,2,tsec)
270  call timab(1506,1,tsec)
271 
272 !===================================
273 !=== Loop on the k-points in IBZ ===
274 !===================================
275  jkg=0
276 
277  if (associated(gs_ham%ph3d_kp)) then
278    nullify (gs_ham%ph3d_kp)
279  end if
280 
281  do jkpt=1,fockbz%mkpt
282 !* nband_k = number of bands at point k_j
283    nband_k=fockbz%nbandocc_bz(jkpt,my_jsppol)
284 !* wtk = weight in BZ of this k point
285    wtk=fockbz%wtk_bz(jkpt) !*sqrt(gs_ham%ucvol)
286 !* jstwfk= how is stored the wavefunction
287    jstwfk=fockbz%istwfk_bz(jkpt)
288 !* npwj= number of plane wave in basis for the wavefunction
289    npwj=fockbz%npwarr(jkpt)
290 !* Basis sphere of G vectors
291    if (allocated(fockbz%cgocc)) then
292      gbound_kp => fockbz%gbound_bz(:,:,jkpt)
293      kg_occ => fockbz%kg_bz(:,1+jkg:npwj+jkg)
294    end if
295 !* Load k^prime hamiltonian in the gs_ham datastructure
296 !  Note: ffnl_kp / ph3d_kp / gbound_kp are not used
297 
298    if (associated(gs_ham%ph3d_kp)) then
299      ABI_ALLOCATE(gs_ham%ph3d_kp,(2,npwj,gs_ham%matblk))
300    end if
301 
302    call load_kprime_hamiltonian(gs_ham,kpt_kp=fockbz%kptns_bz(:,jkpt),&
303 &   istwf_kp=jstwfk,npw_kp=npwj,kg_kp=fockbz%kg_bz(:,1+jkg:npwj+jkg))
304 !* Some temporary allocations needed for PAW
305    if (fockcommon%usepaw==1) then
306      ABI_ALLOCATE(vectin_dum,(2,npwj*nspinor))
307      vectin_dum=zero
308      ABI_ALLOCATE(ffnl_kp_dum,(npwj,0,gs_ham%lmnmax,gs_ham%ntypat))
309      call load_kprime_hamiltonian(gs_ham,ffnl_kp=ffnl_kp_dum)
310    end if
311 
312 ! ======================================
313 ! === Calculate the vector q=k_i-k_j ===
314 ! ======================================
315 !* Evaluation of kpoint_j, the considered k-point in reduced coordinates
316 !     kpoint_j(:)=fockbz%kptns_bz(:,jkpt)
317 !* the vector qvec is expressed in reduced coordinates.
318 !     qvec(:)=kpoint_i(:)-kpoint_j(:)
319    qvec_j(:)=gs_ham%kpt_k(:)-fockbz%kptns_bz(:,jkpt)
320    qeq0=(qvec_j(1)**2+qvec_j(2)**2+qvec_j(3)**2<1.d-15)
321    call bare_vqg(qvec_j,fockcommon%gsqcut,gs_ham%gmet,fockcommon%usepaw,fockcommon%hyb_mixing,&
322 &   fockcommon%hyb_mixing_sr,fockcommon%hyb_range_fock,nfftf,fockbz%nkpt_bz,ngfftf,gs_ham%ucvol,vqg)
323 
324 
325 
326 ! =================================================
327 ! === Loop on the band indices jband of cgocc_k ===
328 ! =================================================
329 
330    do jband=1,nband_k
331 
332 !*   occ = occupancy of jband at this k point
333      occ=fockbz%occ_bz(jband+bdtot_jindex,my_jsppol)
334      if(occ<tol8) cycle
335 
336 ! ==============================================
337 ! === Get cwaveocc_r in real space using FFT ===
338 ! ==============================================
339      if (allocated(fockbz%cwaveocc_bz)) then
340        cwaveocc_r => fockbz%cwaveocc_bz(:,:,:,:,jband+jbg,my_jsppol)
341      else
342        ABI_ALLOCATE(cwaveocc_r,(2,n4f,n5f,n6f))
343        cwaveocc_r=zero
344        call fourwf(1,rhodum0,fockbz%cgocc(:,1+jcg+npwj*(jband-1):jcg+jband*npwj,my_jsppol),rhodum,cwaveocc_r, &
345 &       gbound_kp,gbound_kp,jstwfk,kg_occ,kg_occ,mgfftf,mpi_enreg,ndat1,ngfftf,&
346 &       npwj,1,n4f,n5f,n6f,tim_fourwf0,mpi_enreg%paral_kgb,0,weight1,weight1,use_gpu_cuda=gs_ham%use_gpu_cuda)
347        cwaveocc_r=cwaveocc_r*invucvol
348      end if
349 
350 ! ================================================
351 ! === Get the overlap density matrix rhor_munu ===
352 ! ================================================
353 !* Calculate the overlap density matrix in real space = conj(cwaveocc_r)*cwavef_r
354 !* rhor_munu will contain the overlap density matrix.
355 ! vfock=-int{conj(cwaveocc_r)*cwavef_r*dr'/|r-r'|}
356 
357      call timab(1508,1,tsec)
358      ind=0
359      do i3=1,n3f
360        do i2=1,n2f
361          do i1=1,n1f
362            ind=ind+1
363            recwf  =cwavef_r(1,i1,i2,i3)   ; imcwf  =cwavef_r(2,i1,i2,i3)
364            recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3)
365            rhor_munu(1,ind)= recwocc*recwf+imcwocc*imcwf
366            rhor_munu(2,ind)= recwocc*imcwf-imcwocc*recwf
367          end do ! i1
368        end do ! i2
369      end do ! i3
370      call timab(1508,2,tsec)
371 
372 ! =======================================================
373 ! === Add compensation charge density in the PAW case ===
374 ! === Get the overlap density matrix rhor_munu        ===
375 ! =======================================================
376      call timab(1509,1,tsec)
377      if (fockcommon%usepaw==1) then
378        iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+jband
379 
380        ABI_ALLOCATE(rho12,(2,nfftf,nspinor**2))
381 
382        cwaveocc_prj=>fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1)
383 
384        call pawmknhat_psipsi(cwaveprj,cwaveocc_prj,ider,izero,natom,natom,nfftf,ngfftf,&
385 &       nhat12_grdim,nspinor,fockcommon%ntypat,fockbz%pawang,fockcommon%pawfgrtab,grnhat12,rho12,&
386 &       fockcommon%pawtab,gprimd=gs_ham%gprimd,grnhat_12=grnhat_12,qphon=qvec_j,xred=gs_ham%xred,atindx=gs_ham%atindx)
387 
388        rhor_munu(1,:)=rhor_munu(1,:)+rho12(1,:,nspinor)
389        rhor_munu(2,:)=rhor_munu(2,:)-rho12(2,:,nspinor)
390      end if
391 
392     !Perform an FFT using fourwf to get rhog_munu = FFT^-1(rhor_munu)
393      call fourdp(cplex_fock,rhog_munu,rhor_munu,-1,mpi_enreg,nfftf,ngfftf,mpi_enreg%paral_kgb,tim_fourdp0)
394      call timab(1509,2,tsec)
395 
396      if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
397        call strfock(gs_ham%gprimd,fockcommon%gsqcut,fockstr,fockcommon%hyb_mixing,fockcommon%hyb_mixing_sr,&
398 &       fockcommon%hyb_range_fock,mpi_enreg,nfftf,ngfftf,fockbz%nkpt_bz,rhog_munu,gs_ham%ucvol,qvec_j)
399        fockcommon%stress_ikpt(:,fockcommon%ieigen)=fockcommon%stress_ikpt(:,fockcommon%ieigen)+fockstr(:)*occ*wtk
400        if (fockcommon%usepaw==0.and.(.not.need_ghc)) then
401          if (allocated(fockbz%cgocc)) then
402            ABI_DEALLOCATE(cwaveocc_r)
403          end if
404          cycle
405        end if
406      end if
407 
408 ! ===================================================
409 ! === Calculate the local potential vfockloc_munu ===
410 ! ===================================================
411 !* Apply the Poisson solver to "rhog_munu" while taking into account the effect of the vector "qvec"
412 !* This is precisely what is done in the subroutine hartre, with option cplex=2.
413 !* vfock will contain the local Fock potential, the result of hartre routine.
414 !* vfock = FFT( rhog_munu/|g+qvec|^2 )
415      call timab(1510,1,tsec)
416 #if 0
417 
418      call hartre(cplex_fock,fockcommon%gsqcut,fockcommon%usepaw,mpi_enreg,nfftf,ngfftf,&
419 &     mpi_enreg%paral_kgb,rhog_munu,rprimd,vfock,divgq0=fock%divgq0,qpt=qvec_j)
420 
421 #else
422      do ifft=1,nfftf
423        rhog_munu(1,ifft) = rhog_munu(1,ifft) * vqg(ifft)
424        rhog_munu(2,ifft) = rhog_munu(2,ifft) * vqg(ifft)
425      end do
426 
427 
428      call fourdp(cplex_fock,rhog_munu,vfock,+1,mpi_enreg,nfftf,ngfftf,mpi_enreg%paral_kgb,tim_fourdp0)
429 
430 #endif
431      call timab(1510,2,tsec)
432 
433 !===============================================================
434 !======== Calculate Dij_Fock_hat contribution in case of PAW ===
435 !===============================================================
436 
437      if (fockcommon%usepaw==1) then
438        qphon=qvec_j;nfftotf=product(ngfftf(1:3))
439        cplex_dij=1;ndij=nspden_fock
440        ABI_ALLOCATE(dijhat,(cplex_dij*gs_ham%dimekb1,natom,ndij,cplex_fock))
441        dijhat=zero
442        do iatom=1,natom
443          ipert=iatom
444          itypat=gs_ham%typat(iatom)
445          lmn2_size=fockcommon%pawtab(itypat)%lmn2_size
446          ABI_ALLOCATE(dijhat_tmp,(cplex_fock*cplex_dij*lmn2_size,ndij))
447          dijhat_tmp=zero
448          call pawdijhat(cplex_fock,cplex_dij,dijhat_tmp,gs_ham%gprimd,iatom,ipert,&
449 &         natom,ndij,nfftf,nfftotf,nspden_fock,nspden_fock,fockbz%pawang,fockcommon%pawfgrtab(iatom),&
450 &         fockcommon%pawtab(itypat),vfock,qphon,gs_ham%ucvol,gs_ham%xred)
451          do ii=1,cplex_fock
452            ind=(ii-1)*lmn2_size
453            dijhat(1:cplex_dij*lmn2_size,iatom,:,ii)=dijhat_tmp(ind+1:ind+cplex_dij*lmn2_size,:)
454          end do
455          ABI_DEALLOCATE(dijhat_tmp)
456        end do
457        signs=2; cpopt=2;idir=0; paw_opt=1;nnlout=1;tim_nonlop=1
458        if(need_ghc) then
459          choice=1
460          call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,&
461 &         ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,gvnlc,enl=dijhat,&
462 &         select_k=K_H_KPRIME)
463          ghc2=ghc2-gvnlc*occ*wtk
464        end if
465 
466 ! Forces calculation
467 
468        if (fockcommon%optfor.and.(fockcommon%ieigen/=0)) then
469          choice=2; dotr=zero;doti=zero;cpopt=4
470 
471          do iatom=1,natom
472            do idir=1,3
473              call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,&
474 &             ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,&
475 &             forout,enl=dijhat,iatom_only=iatom,&
476 &             select_k=K_H_KPRIME)
477              call dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,forout,mpi_enreg%me_g0,mpi_enreg%comm_fft)
478              for1(idir)=zero
479              do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd
480                ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft)
481                for1(idir)=for1(idir)+vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-&
482 &               vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom)
483              end do
484            end do
485            do idir=1,3
486              for12(idir)=rprimd(1,idir)*for1(1)+rprimd(2,idir)*for1(2)+rprimd(3,idir)*for1(3)
487              forikpt(idir,iatom)=forikpt(idir,iatom)-(for12(idir)*gs_ham%ucvol/nfftf+dotr(idir))*occ*wtk
488            end do
489          end do
490        end if
491 
492 ! Stresses calculation
493        if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
494          signs=2;choice=3;cpopt=4
495 
496        ! first contribution
497          dotr=zero
498          do idir=1,6
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 dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,strout,mpi_enreg%me_g0,mpi_enreg%comm_fft)
503            fockcommon%stress_ikpt(idir,fockcommon%ieigen)=fockcommon%stress_ikpt(idir,fockcommon%ieigen)-&
504 &           dotr(idir)*occ*wtk/gs_ham%ucvol
505          end do
506        ! second contribution
507          str=zero
508          do iatom=1,natom
509            do idir=1,3
510              do idir1=1,3
511                do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd
512                  ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft)
513                  str(idir,idir1)=str(idir,idir1)+(vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-&
514 &                 vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom))*fockcommon%pawfgrtab(iatom)%rfgd(idir1,ifft)
515 
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_DEALLOCATE(dijhat)
545        ABI_DEALLOCATE(rho12)
546      end if !end PAW
547 
548 ! =============================================================
549 ! === Apply the local potential vfockloc_munu to cwaveocc_r ===
550 ! =============================================================
551      call timab(1507,1,tsec)
552      ind=0
553      do i3=1,ngfftf(3)
554        do i2=1,ngfftf(2)
555          do i1=1,ngfftf(1)
556            ind=ind+1
557 !          ind=i1+ngfftf(1)*(i2-1+ngfftf(2)*(i3-1))
558            revloc=vfock(2*ind-1) ; imvloc=vfock(2*ind)
559            recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3)
560            vlocpsi_r(2*ind-1)=vlocpsi_r(2*ind-1)-(revloc*recwocc-imvloc*imcwocc)*occ*wtk
561            vlocpsi_r(2*ind  )=vlocpsi_r(2*ind  )-(revloc*imcwocc+imvloc*recwocc)*occ*wtk
562          end do
563        end do
564      end do
565      call timab(1507,2,tsec)
566      if (allocated(fockbz%cgocc)) then
567        ABI_DEALLOCATE(cwaveocc_r)
568      end if
569 
570    end do ! jband
571 
572 ! ================================================
573 ! === End : update of shifts and deallocations ===
574 ! ================================================
575 !* Update of the shifts to be applied (reminder : mkmem is not 0, nspinor=1)
576    jcg=jcg+npwj*nband_k
577    jbg=jbg+nband_k
578    bdtot_jindex=bdtot_jindex+nband_k
579    jkg=jkg+npwj
580    if (fockcommon%usepaw==1) then
581      ABI_DEALLOCATE(vectin_dum)
582      ABI_DEALLOCATE(ffnl_kp_dum)
583    end if
584    if (associated(gs_ham%ph3d_kp)) then
585      ABI_DEALLOCATE(gs_ham%ph3d_kp)
586    end if
587  end do ! jkpt
588 
589  if (fockcommon%usepaw==1) then
590    if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
591      call xmpi_sum(forikpt,mpi_enreg%comm_hf,ier)
592      do iatom=1,natom !Loop over atom
593        ia=gs_ham%atindx(iatom)
594        fockcommon%forces_ikpt(:,ia,fockcommon%ieigen)=forikpt(:,iatom)
595      end do
596    end if
597  end if
598  if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
599    call xmpi_sum(fockcommon%stress_ikpt,mpi_enreg%comm_hf,ier)
600  end if
601 
602  if (.not.need_ghc) then
603 
604 ! ===============================
605 ! === Deallocate local arrays ===
606 ! ===============================
607    ABI_DEALLOCATE(cwavef_r)
608    ABI_DEALLOCATE(ghc1)
609    ABI_DEALLOCATE(ghc2)
610    ABI_DEALLOCATE(rhor_munu)
611    ABI_DEALLOCATE(rhog_munu)
612    ABI_DEALLOCATE(vlocpsi_r)
613    ABI_DEALLOCATE(dummytab)
614    ABI_DEALLOCATE(vfock)
615    ABI_DEALLOCATE(vqg)
616    if (fockcommon%usepaw==1) then
617      ABI_DEALLOCATE(gvnlc)
618      ABI_DEALLOCATE(grnhat12)
619      if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
620        ABI_DEALLOCATE(forikpt)
621        ABI_DEALLOCATE(forout)
622      end if
623      if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
624        ABI_DEALLOCATE(strout)
625      end if
626      ABI_DEALLOCATE(grnhat_12)
627    end if
628    if(fockcommon%usepaw==1.or.fockcommon%optstr) then
629      ABI_DEALLOCATE(gboundf)
630    end if
631 !*Restore gs_ham datastructure
632 
633    if (associated(gs_ham%ph3d_kp)) then
634      ABI_ALLOCATE(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk))
635    end if
636    call load_kprime_hamiltonian(gs_ham,kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,&
637 &   npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k)
638 
639 !   if (fockcommon%ieigen/=0) fockcommon%ieigen=0
640    return
641  end if
642 
643 
644  call timab(1506,2,tsec)
645  call timab(1511,1,tsec)
646 
647 !*Restore gs_ham datastructure
648 
649  if (associated(gs_ham%ph3d_kp)) then
650    ABI_ALLOCATE(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk))
651  end if
652  call load_kprime_hamiltonian(gs_ham,kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,&
653 & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k)
654 
655 !* Perform an FFT using fourwf to get ghc1 = FFT^-1(vlocpsi_r)
656  ABI_ALLOCATE(psilocal,(cplex_fock*n4f,n5f,n6f))
657  call fftpac(1,mpi_enreg,nspden_fock,cplex_fock*n1f,n2f,n3f,cplex_fock*n4f,n5f,n6f,ngfft,vlocpsi_r,psilocal,2)
658 
659  call fourwf(0,rhodum0,rhodum,ghc1,psilocal,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,&
660 & mgfftf,mpi_enreg,ndat1,ngfftf,1,npw,n4f,n5f,n6f,3,mpi_enreg%paral_kgb,tim_fourwf0,weight1,weight1,&
661 & use_gpu_cuda=gs_ham%use_gpu_cuda)
662  ABI_DEALLOCATE(psilocal)
663 
664  ghc1=ghc1*sqrt(gs_ham%ucvol)+ghc2
665 
666 !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc
667  ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_hf + ghc1(:,:)
668 
669  call xmpi_sum(ghc,mpi_enreg%comm_hf,ier)
670 
671  call timab(1511,2,tsec)
672 
673 
674 ! ===============================
675 ! === Deallocate local PAW arrays ===
676 ! ===============================
677 
678  if (fockcommon%usepaw==1) then
679    ABI_DEALLOCATE(gvnlc)
680    ABI_DEALLOCATE(grnhat12)
681    if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
682      ABI_DEALLOCATE(forikpt)
683      ABI_DEALLOCATE(forout)
684    end if
685    if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
686      ABI_DEALLOCATE(strout)
687    end if
688    ABI_DEALLOCATE(grnhat_12)
689  end if
690  if(fockcommon%usepaw==1.or.fockcommon%optstr) then
691    ABI_DEALLOCATE(gboundf)
692  end if
693 
694 ! ============================================
695 ! === Calculate the contribution to energy ===
696 ! ============================================
697 !* Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction
698 !* at each convergence step. (consistent definition with the defintion of hartree energy)
699  if (fockcommon%ieigen/=0) then
700    eigen=zero
701 !* Dot product of cwavef and ghc
702 !* inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering
703    if(gs_ham%istwf_k==2) then
704      eigen=half*cwavef(1,1)*ghc1(1,1)
705    else
706      eigen=cwavef(1,1)*ghc1(1,1)+cwavef(2,1)*ghc1(2,1)
707    end if
708    do ipw=2,npw
709      eigen=eigen+cwavef(1,ipw)*ghc1(1,ipw)+cwavef(2,ipw)*ghc1(2,ipw)
710    end do
711    if(gs_ham%istwf_k>=2) eigen=two*eigen
712    call xmpi_sum(eigen,mpi_enreg%comm_hf,ier)
713    fockcommon%eigen_ikpt(fockcommon%ieigen)= eigen
714    if(fockcommon%use_ACE==0) fockcommon%ieigen = 0
715  end if
716 
717 ! ===============================
718 ! === Deallocate local arrays ===
719 ! ===============================
720  ABI_DEALLOCATE(cwavef_r)
721  ABI_DEALLOCATE(ghc1)
722  ABI_DEALLOCATE(ghc2)
723  ABI_DEALLOCATE(rhor_munu)
724  ABI_DEALLOCATE(rhog_munu)
725  ABI_DEALLOCATE(vlocpsi_r)
726  ABI_DEALLOCATE(dummytab)
727  ABI_DEALLOCATE(vfock)
728  ABI_DEALLOCATE(vqg)
729 
730  call timab(1504,2,tsec)
731 
732 end subroutine fock_getghc

ABINIT/m_fock_getghc [ Modules ]

[ Top ] [ Modules ]

NAME

FUNCTION

COPYRIGHT

  Copyright (C) 2013-2018 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 .

PARENTS

CHILDREN

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_fock_getghc
26 
27  use defs_basis
28  use defs_abitypes
29  use m_abicore
30  use m_errors
31  use m_xmpi
32  use m_fock
33  use m_pawcprj
34  !use m_cgtools
35 
36  use defs_datatypes, only : pseudopotential_type
37  use m_time,         only : timab
38  use m_symtk,        only : matr3inv
39  use m_cgtools,      only : dotprod_g
40  use m_kg,           only : mkkpg
41  use m_fftcore,      only : sphereboundary
42  use m_fft,          only : fftpac, fourwf, fourdp
43  use m_hamiltonian,  only : gs_hamiltonian_type,load_kprime_hamiltonian,K_H_KPRIME,load_k_hamiltonian, &
44                             init_hamiltonian, destroy_hamiltonian, load_spin_hamiltonian
45  use m_pawdij,       only : pawdijhat
46  use m_pawrhoij,     only : pawrhoij_type, pawrhoij_free, pawrhoij_alloc
47  use m_paw_nhat,     only : pawmknhat_psipsi
48  use m_spacepar,     only : hartre
49  use m_nonlop,       only : nonlop
50  use m_bandfft_kpt,      only : bandfft_kpt, bandfft_kpt_type, bandfft_kpt_savetabs,bandfft_kpt_restoretabs, &
51                                 prep_bandfft_tabs
52  use m_pawtab,           only : pawtab_type
53  use m_paw_ij,           only : paw_ij_type
54  use m_mkffnl,           only : mkffnl
55  use m_mpinfo,           only : proc_distrb_cycle
56 
57  implicit none
58 
59  private