TABLE OF CONTENTS


ABINIT/pawnhatfr [ Functions ]

[ Top ] [ Functions ]

NAME

 pawnhatfr

FUNCTION

 PAW: Compute frozen part of charge compensation density nhat
      nhatfr(r)=Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
      Depends on q wave vector but not on first-order wave-function.

COPYRIGHT

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

INPUTS

  ider=0: computes frozen part of compensation density
       1: computes frozen part of compensation density and cartesian gradients
  idir=direction of atomic displacement (in case of phonons perturb.)
  ipert=nindex of perturbation
  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
  natom=total number of atoms in cell
  nspden=number of spin-density components
  ntypat=number of types of atoms
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrhoij(my_natom) <type(pawrhoij_type)>= Ground-State paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  rprimd(3,3)=dimensional primitive translations for real space

OUTPUT

  pawfgrtab(iatom)%nhatfr(nfgd,nspden)
                  frozen part of charge compensation density (inside PAW spheres)
                  =Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
  === If ider==1
  pawfgrtab(iatom)%nhatfrgr(3,nfgd,nspden)
                  gradients of frozen part of charge compensation density (inside PAW spheres)
                  =Sum_ij,lm[rhoij_ij.q_ij^l . d/dr((g_l(r).Y_lm(r))^(1))]

PARENTS

      dfpt_nstpaw,dfpt_scfcv,pawmknhat

CHILDREN

      free_my_atmtab,get_my_atmtab,pawgylm

SOURCE

 53 #if defined HAVE_CONFIG_H
 54 #include "config.h"
 55 #endif
 56 
 57 #include "abi_common.h"
 58 
 59 subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,&
 60 &                    pawang,pawfgrtab,pawrhoij,pawtab,rprimd, &
 61 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
 62 
 63  use defs_basis
 64  use m_profiling_abi
 65  use m_errors
 66  use m_xmpi, only : xmpi_comm_self
 67 
 68  use m_pawang,       only : pawang_type
 69  use m_pawtab,       only : pawtab_type
 70  use m_pawfgrtab,    only : pawfgrtab_type
 71  use m_pawrhoij,     only : pawrhoij_type
 72  use m_paw_finegrid, only : pawgylm
 73  use m_paral_atom,   only : get_my_atmtab, free_my_atmtab
 74 
 75 !This section has been created automatically by the script Abilint (TD).
 76 !Do not modify the following lines by hand.
 77 #undef ABI_FUNC
 78 #define ABI_FUNC 'pawnhatfr'
 79 !End of the abilint section
 80 
 81  implicit none
 82 
 83 !Arguments ------------------------------------
 84 !scalars
 85  integer,intent(in) :: ider,idir,ipert,my_natom,natom,nspden,ntypat
 86  integer,optional,intent(in) :: comm_atom
 87  type(pawang_type),intent(in) :: pawang
 88 !arrays
 89  integer,optional,target,intent(in) :: mpi_atmtab(:)
 90  real(dp),intent(in) :: rprimd(3,3)
 91  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
 92  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
 93  type(pawtab_type),intent(in) :: pawtab(ntypat)
 94 
 95 !Local variables-------------------------------
 96 !scalars
 97  integer :: iatom,iatom_tot,ic,ils,ilslm,irhoij,isel,ispden,istr,itypat,jrhoij
 98  integer :: klm,klmn,lm_size,lmn2_size,lm0,lmax,lmin,mua,mub,mm,mu,my_comm_atom,nfgd,nu,optgr0,optgr1,optgr2
 99  logical :: my_atmtab_allocated,my_pert,paral_atom
100  real(dp) :: contrib,ro
101 !arrays
102  integer,parameter :: voigt(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
103  integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
104  integer,pointer :: my_atmtab(:)
105  real(dp),allocatable :: nhatfr_tmp(:,:),nhatfrgr_tmp(:,:,:)
106 
107 ! *************************************************************************
108 
109  DBG_ENTER("COLL")
110 
111 !Only relevant for atomic displacement and strain perturbation
112  if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4) return
113 
114 !Set up parallelism over atoms
115  paral_atom=(present(comm_atom).and.(my_natom/=natom))
116  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
117  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
118  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
119 
120 !Compatibility tests
121  if (my_natom>0) then
122    if ((pawfgrtab(1)%gylm_allocated==0.or.pawfgrtab(1)%gylmgr_allocated==0).and. &
123 &   pawfgrtab(1)%rfgd_allocated==0) then
124      MSG_BUG('  pawfgrtab()%rfgd array must be allocated  !')
125    end if
126  end if
127 
128  my_pert = (ipert<=natom).or.ipert==natom+3.or.ipert==natom+4
129 
130 !Get correct index of strain pertubation
131  if (ipert==natom+3) istr = idir
132  if (ipert==natom+4) istr = idir + 3
133  
134 !Loops over  atoms
135  do iatom=1,my_natom
136    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
137 
138 !  Eventually allocate frozen nhat points
139    if (my_pert) then
140      if (pawfgrtab(iatom)%nhatfr_allocated==0) then
141        if (allocated(pawfgrtab(iatom)%nhatfr))  then
142          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
143        end if
144        ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(pawfgrtab(iatom)%nfgd,nspden))
145        pawfgrtab(iatom)%nhatfr_allocated=1
146      end if
147      if (ider==1.and.pawfgrtab(iatom)%nhatfrgr_allocated==0) then
148        if (allocated(pawfgrtab(iatom)%nhatfrgr))  then
149          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr)
150        end if
151        ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(3,pawfgrtab(iatom)%nfgd,nspden))
152        pawfgrtab(iatom)%nhatfrgr_allocated=1
153      end if
154    end if
155 
156 !  Select if frozen part of nhat exists for the current perturbation
157    if ((.not.my_pert).or.(pawfgrtab(iatom)%nhatfr_allocated==0)) cycle
158 
159 !  Some atom-dependent quantities
160    itypat=pawfgrtab(iatom)%itypat
161    lm_size=pawfgrtab(iatom)%l_size**2
162    lmn2_size=pawtab(itypat)%lmn2_size
163 
164 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
165    nfgd=pawfgrtab(iatom)%nfgd
166    if ((pawfgrtab(iatom)%gylmgr_allocated==0).or. &
167 &   (pawfgrtab(iatom)%gylmgr2_allocated==0.and.ider==1)) then
168      optgr0=0;optgr1=0;optgr2=0
169      if(ipert==natom+3.or.ipert==natom+4)then
170        if (pawfgrtab(iatom)%gylm_allocated==0) then
171          if (allocated(pawfgrtab(iatom)%gylm))  then
172            ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
173          end if
174          ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
175          pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
176        end if
177      end if
178      if (pawfgrtab(iatom)%gylmgr_allocated==0) then
179        if (allocated(pawfgrtab(iatom)%gylmgr))  then
180          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
181        end if
182        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size))
183        pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
184      end if
185      if (ider==1.and.pawfgrtab(iatom)%gylmgr2_allocated==0) then
186        if (allocated(pawfgrtab(iatom)%gylmgr2))  then
187          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
188        end if
189        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size))
190        pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1
191      end if
192      if (optgr0+optgr1+optgr2>0) then
193        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
194 &       lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd)
195      end if
196    end if
197 
198 
199 !  ============ Phonons ====================================
200    if (ipert<=natom) then
201 
202 !    Loop over spin components
203      do ispden=1,nspden
204 
205        ABI_ALLOCATE(nhatfr_tmp,(3,nfgd))
206        nhatfr_tmp=zero
207        if (ider==1) then
208          ABI_ALLOCATE(nhatfrgr_tmp,(3,nfgd,3))
209          nhatfrgr_tmp=zero
210        end if
211        
212        jrhoij=1
213        do irhoij=1,pawrhoij(iatom)%nrhoijsel
214          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
215          klm =pawtab(itypat)%indklmn(1,klmn)
216          lmin=pawtab(itypat)%indklmn(3,klmn)
217          lmax=pawtab(itypat)%indklmn(4,klmn)
218          if (nspden/=2) then
219            ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
220          else
221            if (ispden==1) then
222              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
223            else if (ispden==2) then
224              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
225            end if
226          end if
227          ro=pawtab(itypat)%dltij(klmn)*ro
228          do ils=lmin,lmax,2
229            lm0=ils**2+ils+1
230            do mm=-ils,ils
231              ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
232              if (isel>0) then
233                do ic=1,nfgd
234                  do mu=1,3
235                    contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn)&
236 &                   *pawfgrtab(iatom)%gylmgr(mu,ic,ilslm)
237                    nhatfr_tmp(mu,ic)=nhatfr_tmp(mu,ic)+contrib
238                  end do
239                end do
240                if (ider==1) then
241                  do ic=1,nfgd
242                    do nu=1,3
243                      do mu=1,3
244                        contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
245 &                       *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
246                        nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
247                      end do
248                    end do
249                  end do
250                end if
251              end if
252            end do
253          end do
254          jrhoij=jrhoij+pawrhoij(iatom)%cplex
255        end do
256 
257 !      Convert from cartesian to reduced coordinates
258        do ic=1,nfgd
259          pawfgrtab(iatom)%nhatfr(ic,ispden)= &
260 &         rprimd(1,idir)*nhatfr_tmp(1,ic) &
261 &         +rprimd(2,idir)*nhatfr_tmp(2,ic) &
262 &         +rprimd(3,idir)*nhatfr_tmp(3,ic)
263        end do
264        if (ider==1) then
265          do nu=1,3
266            do ic=1,nfgd
267              pawfgrtab(iatom)%nhatfrgr(nu,ic,ispden)= &
268 &             rprimd(1,idir)*nhatfrgr_tmp(1,ic,nu) &
269 &             +rprimd(2,idir)*nhatfrgr_tmp(2,ic,nu) &
270 &             +rprimd(3,idir)*nhatfrgr_tmp(3,ic,nu)
271            end do
272          end do
273        end if
274        ABI_DEALLOCATE(nhatfr_tmp)
275        if (ider==1) then
276          ABI_DEALLOCATE(nhatfrgr_tmp)
277        end if
278 !      End loop over spin components
279      end do ! ispden
280 
281 
282 !  ============ Elastic tensor ===============================
283    else if (ipert==natom+3.or.ipert==natom+4) then
284 !    Loop over spin components
285      pawfgrtab(iatom)%nhatfr(:,:) = zero
286      do ispden=1,nspden
287        jrhoij=1
288        do irhoij=1,pawrhoij(iatom)%nrhoijsel
289          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
290          klm =pawtab(itypat)%indklmn(1,klmn)
291          lmin=pawtab(itypat)%indklmn(3,klmn)
292          lmax=pawtab(itypat)%indklmn(4,klmn)
293          if (nspden/=2) then
294            ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
295          else
296            if (ispden==1) then
297              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
298            else if (ispden==2) then
299              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
300            end if
301          end if
302          ro=pawtab(itypat)%dltij(klmn)*ro
303          do ils=lmin,lmax,2
304            lm0=ils**2+ils+1
305            do mm=-ils,ils
306              ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
307              if (isel>0) then
308 !              Sum{[Q_ij_q^LM^(1)]}
309                do ic=1,nfgd
310                  mua=alpha(istr);mub=beta(istr)
311                  pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
312 &                 ro*pawtab(itypat)%qijl(ilslm,klmn)*half*(&
313 &                 pawfgrtab(iatom)%gylmgr(mua,ic,ilslm)*pawfgrtab(iatom)%rfgd(mub,ic)&
314 &                 +pawfgrtab(iatom)%gylmgr(mub,ic,ilslm)*pawfgrtab(iatom)%rfgd(mua,ic))
315                end do
316 !              Add volume contribution
317                if(istr<=3)then
318                  do ic=1,nfgd
319                    pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
320 &                   ro*pawtab(itypat)%qijl(ilslm,klmn)*pawfgrtab(iatom)%gylm(ic,ilslm)
321                  end do
322                end if
323                if (ider==1) then
324                  MSG_ERROR("nhatgr not implemented for strain perturbationxs")
325 !                 do ic=1,nfgd
326 !                   do nu=1,6
327 !                     do mu=1,6
328 !                       contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
329 !&                       *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
330 !                       nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
331 !                     end do
332 !                   end do
333 !                 end do
334                end if
335              end if
336            end do
337          end do
338          jrhoij=jrhoij+pawrhoij(iatom)%cplex
339        end do
340      end do ! ispden
341    end if
342 
343 !  Eventually free temporary space for g_l(r).Y_lm(r) gradients and exp(-i.q.r)
344    if (pawfgrtab(iatom)%gylmgr_allocated==2) then
345      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
346      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0))
347      pawfgrtab(iatom)%gylmgr_allocated=0
348    end if
349    if (pawfgrtab(iatom)%gylmgr2_allocated==2) then
350      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
351      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(0,0,0))
352      pawfgrtab(iatom)%gylmgr2_allocated=0
353    end if
354 
355 !  End loop on atoms
356  end do
357 
358 
359 !Destroy atom table used for parallelism
360  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
361 
362  DBG_EXIT("COLL")
363 
364 end subroutine pawnhatfr