TABLE OF CONTENTS


ABINIT/paw_mknewh0 [ Functions ]

[ Top ] [ Functions ]

NAME

 paw_mknewh0

FUNCTION

 Calculates the new bare PAW Hamiltonian in the case of quasi-particle self-consistent GW calculations.

COPYRIGHT

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

INPUTS

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  nsppol=1 for unpolarized, 2 for spin-polarized
  nspden=number of spin-density components
  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid
  pawspnorb=flag: 1 if spin-orbit coupling is activated
  pawprtvol=control print volume and debugging output for PAW
  Cryst<crystal_t>=Info on unit cell and its symmetries
  Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
  Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
  Pawang<type(pawang_type)>=paw angular mesh and related data
  Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  vxc(nfftf,nspden)=exchange-correlation potential
  vxc_val(nfftf,nspden)=valence only exchange-correlation potential
  vtrial(nfftf,nspden)=potential (Hartree+XC+loc)

SIDE EFFECTS

  Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
     At output: new value for Paw_ij()%dij

PARENTS

      calc_vhxc_me

CHILDREN

      free_my_atmtab,get_my_atmtab,pawgylm,symdij,symdij_all,wrtout

SOURCE

 44 #if defined HAVE_CONFIG_H
 45 #include "config.h"
 46 #endif
 47 
 48 #include "abi_common.h"
 49 
 50 subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,&
 51 &          Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,&
 52 &          mpi_atmtab,comm_atom) ! optional arguments (parallelism)
 53 
 54  use defs_basis
 55  use m_profiling_abi
 56  use m_errors
 57  use m_xmpi, only : xmpi_comm_self
 58 
 59  use m_crystal,      only : crystal_t
 60  use m_pawang,       only : pawang_type
 61  use m_pawtab,       only : pawtab_type
 62  use m_paw_an,       only : paw_an_type
 63  use m_paw_ij,       only : paw_ij_type
 64  use m_pawfgrtab,    only : pawfgrtab_type
 65  use m_pawdij,       only : symdij, symdij_all
 66  use m_paw_finegrid, only : pawgylm
 67  use m_paral_atom,   only : get_my_atmtab, free_my_atmtab
 68 
 69 !This section has been created automatically by the script Abilint (TD).
 70 !Do not modify the following lines by hand.
 71 #undef ABI_FUNC
 72 #define ABI_FUNC 'paw_mknewh0'
 73  use interfaces_14_hidewrite
 74 !End of the abilint section
 75 
 76  implicit none
 77 
 78 !Arguments ------------------------------------
 79 !scalars
 80  integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb
 81  integer,optional,intent(in) :: comm_atom
 82 !arrays
 83  integer,optional,target,intent(in) :: mpi_atmtab(:)
 84  real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden)
 85  type(crystal_t),intent(in) :: Cryst
 86  type(Pawang_type),intent(in) :: Pawang
 87  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
 88  type(Paw_an_type),intent(in) :: Paw_an(my_natom)
 89  type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom)
 90  type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom)
 91 
 92 !Local variables-------------------------------
 93 !scalars
 94  integer,parameter :: ipert0=0
 95  integer :: iat,iat_tot,idij,cplex,ndij,option_dij
 96  integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm
 97  integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij
 98  integer :: ils,ilslm,ic,lm0
 99  integer :: nsploop,is2fft
100  real(dp) :: gylm,qijl
101  logical :: ltest,my_atmtab_allocated,paral_atom
102  character(len=500) :: msg
103 !arrays
104  integer, pointer :: indklmn_(:,:)
105  integer,pointer :: my_atmtab(:)
106  real(dp) :: rdum(1),rdum2(1)
107  real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:)
108  real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:)
109 
110 ! *************************************************************************
111 
112  DBG_ENTER("COLL")
113 
114  call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL')
115 
116 !== Set up parallelism over atoms ===
117  paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom))
118  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
119  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
120  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom)
121 
122  if (my_natom>0) then
123 
124 !  === Test if required pointers in paw_ij are allocated ===
125    ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) )
126    ABI_CHECK(ltest,'dijxc or dijxc_val not calculated')
127 
128    ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2)
129    ABI_CHECK(ltest,'dijhat not calculated')
130 
131    ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2)
132    ABI_CHECK(ltest,'dijhartree not calculated')
133 
134    if (ANY(Pawtab(:)%usepawu>0)) then
135      do iat=1,my_natom
136        iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
137        itypat=Cryst%typat(iat_tot)
138        if (Pawtab(itypat)%usepawu>0) then
139          ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2)
140          write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated."
141          ABI_CHECK(ltest,msg)
142        end if
143      end do
144    end if
145 
146    if (pawspnorb>0) then
147      do iat=1,my_natom
148        ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2)
149        write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated."
150        ABI_CHECK(ltest,msg)
151      end do
152    end if
153  end if ! my_natom>0
154 
155 !== Construct the new PAW H0 Hamiltonian ===
156  do iat=1,my_natom
157    iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
158 
159    itypat    = Cryst%typat(iat_tot)
160    lmn_size  = Pawtab(itypat)%lmn_size
161    lmn2_size = Pawtab(itypat)%lmn2_size
162    lm_size   = Paw_an(iat)%lm_size
163    cplex     = Paw_ij(iat)%cplex
164    cplex_dij = Paw_ij(iat)%cplex_dij
165    ndij      = Paw_ij(iat)%ndij
166 
167    ABI_CHECK(cplex==1,'cplex/=1 not implemented')
168    ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented')
169 !  
170 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
171    if (Pawfgrtab(iat)%gylm_allocated==0) then
172      if (allocated(Pawfgrtab(iat)%gylm))  then
173        ABI_DEALLOCATE(Pawfgrtab(iat)%gylm)
174      end if
175      ABI_ALLOCATE(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size))
176      Pawfgrtab(iat)%gylm_allocated=2
177 
178      call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,&
179 &     Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd)
180    end if
181 
182 !  === Calculate LM contribution to dijhmxc_val for this atom ===
183 !  * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not
184 !  contain neither the hat term nor the LM sum of onsite terms (they should cancel each other)
185 !  FIXME change paw_dij,  otherwise I miss tnc in vxc
186 !  * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh ===
187 !  * The following quantities do not depend on ij
188    ABI_ALLOCATE(prod_hloc   ,(lm_size,ndij))
189    ABI_ALLOCATE(prodhxc_core,(lm_size,ndij))
190    prod_hloc   =zero
191    prodhxc_core=zero
192    do idij=1,ndij
193      do ilslm=1,lm_size
194        do ic=1,Pawfgrtab(iat)%nfgd
195          is2fft=Pawfgrtab(iat)%ifftsph(ic)
196          gylm=Pawfgrtab(iat)%gylm(ic,ilslm)
197          prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm
198 !        prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm
199          prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm
200        end do
201      end do
202    end do !idij
203 
204 !  === Assembly the "Hat" contribution for this atom ====
205    ABI_ALLOCATE(dijhl_hat  ,(cplex_dij*lmn2_size,ndij))
206    ABI_ALLOCATE(dijhmxc_val,(cplex_dij*lmn2_size,ndij))
207    dijhl_hat  =zero
208    dijhmxc_val=zero
209    indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size)
210 
211    do idij=1,ndij
212      do klmn=1,lmn2_size
213        klm =indklmn_(1,klmn)
214        lmin=indklmn_(3,klmn)
215        lmax=indklmn_(4,klmn)
216 
217 !      === $\sum_lm q_ij^l prod* for each idij$ ===
218        do ils=lmin,lmax,2
219          lm0=ils**2+ils+1
220          do mm=-ils,ils
221            ilslm=lm0+mm
222            isel=Pawang%gntselect(lm0+mm,klm)
223            if (isel>0) then
224              qijl=Pawtab(itypat)%qijl(ilslm,klmn)
225              dijhl_hat  (klmn,idij)=dijhl_hat  (klmn,idij) +  prod_hloc (ilslm,idij)*qijl
226              dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl
227            end if
228          end do
229        end do
230      end do
231    end do
232 
233    ABI_DEALLOCATE(prod_hloc)
234    ABI_DEALLOCATE(prodhxc_core)
235 
236 !  * Normalization factor due to integration on the FFT mesh
237    dijhl_hat  = dijhl_hat  *Cryst%ucvol/DBLE(nfftf)
238    dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf)
239 
240 !  === Now assembly the bare Hamiltonian ===
241 !  * Loop over density components overwriting %dij
242    nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4
243 
244    do idij=1,nsploop
245      klmn1=1
246 
247      do jlmn=1,lmn_size
248        j0lmn=jlmn*(jlmn-1)/2
249        do ilmn=1,jlmn
250          klmn=j0lmn+ilmn
251 
252 !        The following gives back the input dij.
253 !        since dijxc contains the hat term done on the FFT mesh
254          if (.FALSE.) then
255            Paw_ij(iat)%dij(klmn,idij) =        &
256 &           Pawtab(itypat)%dij0    (klmn)      &
257 &           +Paw_ij(iat)%dijhartree(klmn)      &
258 &           +Paw_ij(iat)%dijxc     (klmn,idij) &
259 &           +dijhl_hat   (klmn,idij)
260 
261          else
262 !          === Make nonlocal part of h0 removing the valence contribution ===
263 !          Remeber that XC contains already the Hat contribution
264            Paw_ij(iat)%dij(klmn,idij) =        &
265 &           Pawtab(itypat)%dij0      (klmn)    &
266 &           +Paw_ij(iat)%dijhartree(klmn)      &
267 &           +Paw_ij(iat)%dijxc     (klmn,idij) &  ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val)
268 &          -Paw_ij(iat)%dijxc_val (klmn,idij) &  ! I suppose that the "hat" term on the FFT mesh in included in both.
269 &          +dijhmxc_val(klmn,idij)               ! Local + Hartree - XC val contribution to the "hat" term.
270 
271 !          Add the U contribution to the 
272 !          if (.FALSE. .and. Pawtab(itypat)%usepawu>0) then
273            if (.TRUE. .and. Pawtab(itypat)%usepawu>0) then
274              Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij)
275            end if
276          end if
277 !        TODO dijso, dijU, vpawx?
278 !        Just to be consistent, update some values.
279 !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij)
280 
281        end do !ilmn
282      end do !jlmn
283    end do !idij
284 
285 !  this is to be consistent?
286 !  deallocate(Paw_ij(iat)%dijvxc_val)
287    ABI_DEALLOCATE(dijhl_hat)
288    ABI_DEALLOCATE(dijhmxc_val)
289  end do !iat
290 
291 !=== Symmetrize total Dij ===
292  option_dij=0 ! For total Dij.
293 #if 0
294  if (paral_atom) then
295    call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
296 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
297 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
298  else
299    call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
300 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
301  end if
302 #else
303  if (paral_atom) then
304    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
305 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
306 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
307  else
308    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
309 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
310  end if
311 #endif
312 
313 !Destroy atom table used for parallelism
314  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
315 
316  DBG_EXIT("COLL")
317 
318 end subroutine paw_mknewh0