TABLE OF CONTENTS


ABINIT/m_dens [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dens

FUNCTION

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group ()
 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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_dens
27 
28  use defs_basis
29  use m_errors
30  use m_profiling_abi
31  use m_splines
32 
33  use m_geometry,         only : xcart2xred, metric
34 
35  implicit none
36 
37  private
38 
39 ! public procedures.
40  public :: dens_hirsh      ! Compute the Hirshfeld charges

m_dens/dens_hirsh [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 dens_hirsh

FUNCTION

 Compute the Hirshfeld charges

INPUTS

  mpoint=Maximum number of points in radial meshes.
  radii(mpoint, ntypat)=Radial meshes for each type
  aeden(mpoint, nytpat)=All-electron densities.
  npoint(ntypat)=The number of the last point with significant density is stored in npoint(itypat)
  minimal_den=Tolerance on the minum value of the density
  grid_den(nrx,nry,nrz)= density on the grid
  natom = number of atoms in the unit cell
  nrx,nry,nrz= number of points in the grid for the three directions
  ntypat=number of types of atoms in unit cell.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  typat(natom)=type of each atom
  xcart(3,natom) = different positions of the atoms in the unit cell
  zion=(ntypat)gives the ionic charge for each type of atom
  prtcharge=1 to write the Hirshfeld charge decomposition

OUTPUT

  hcharge(natom), hden(natom), hweight(natom)= Hirshfeld charges, densities, weights.

PARENTS

      m_cut3d

CHILDREN

      metric,spline,xcart2xred

SOURCE

 81 subroutine dens_hirsh(mpoint,radii,aeden,npoint,minimal_den,grid_den, &
 82   natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,prtcharge,hcharge,hden,hweight)
 83 
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'dens_hirsh'
 89 !End of the abilint section
 90 
 91  implicit none
 92 
 93 !Arguments ------------------------------------
 94 !scalars
 95  integer,intent(in) :: natom,nrx,nry,nrz,ntypat,prtcharge,mpoint
 96  real(dp),intent(in) :: minimal_den
 97 !arrays
 98  integer,intent(in) :: typat(natom),npoint(ntypat)
 99  real(dp),intent(in) :: grid_den(nrx,nry,nrz),rprimd(3,3),zion(ntypat)
100  real(dp),intent(in) :: xcart(3,natom)
101  real(dp),intent(in) :: radii(mpoint,ntypat),aeden(mpoint,ntypat)
102  real(dp),intent(out) :: hcharge(natom),hden(natom),hweight(natom)
103 
104 !Local variables -------------------------
105 !scalars
106  integer :: i1,i2,i3,iatom,icell,igrid,ii,inmax,inmin,istep,itypat
107  integer :: k1,k2,k3,mcells,nfftot,ngoodpoints,npt
108  real(dp) :: aa,bb,coeff1,coeff2,coeff3,den,factor,h_inv,hh,maxrad
109  real(dp) :: rr,rr2,total_charge,total_weight,total_zion,ucvol
110  real(dp) :: yp1,ypn
111  !character(len=500) :: msg
112 !arrays
113  integer :: highest(3),lowest(3)
114  integer,allocatable :: ncells(:)
115  real(dp) :: coordat(3),coord23_1,coord23_2,coord23_3,diff1,diff2,diff3,gmet(3,3),gprimd(3,3),rmet(3,3)
116  real(dp) :: vperp(3),width(3)
117  real(dp),allocatable :: coord1(:,:),local_den(:,:,:,:)
118  real(dp),allocatable :: step(:,:),sum_den(:,:,:),work(:)
119  real(dp),allocatable :: xcartcells(:,:,:),xred(:,:),yder2(:)
120 
121 ! *********************************************************************
122 
123 !1. Read the 1D all-electron atomic files
124 !Store the radii in radii(:,itypat), and the all-electron
125 !densities in aeden(:,itypat). The number of the last
126 !point with significant density is stored in npoint(itypat)
127 
128 !2. Compute the list of atoms that are sufficiently close to the cell
129 
130  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
131  nfftot=nrx*nry*nrz
132 !DEBUG
133 !do k3=1,nrz
134 !do k2=1,nry
135 !do k1=1,nrx
136 !total_charge=total_charge+grid_den(k1,k2,k3)
137 !end do
138 !end do
139 !end do
140 !write(std_out,*)' total_charge=',total_charge*ucvol/dble(nfftot)
141 !ENDDEBUG
142 
143  ABI_ALLOCATE(xred,(3,natom))
144  call xcart2xred(natom,rprimd,xcart,xred)
145 
146 !Compute the widths of the cell
147 !First width : perpendicular vector length
148  vperp(:)=rprimd(:,1)-rprimd(:,2)*rmet(1,2)/rmet(2,2) -rprimd(:,3)*rmet(1,3)/rmet(3,3)
149  width(1)=sqrt(dot_product(vperp,vperp))
150 !Second width
151  vperp(:)=rprimd(:,2)-rprimd(:,1)*rmet(2,1)/rmet(1,1) -rprimd(:,3)*rmet(2,3)/rmet(3,3)
152  width(2)=sqrt(dot_product(vperp,vperp))
153 !Third width
154  vperp(:)=rprimd(:,3)-rprimd(:,1)*rmet(3,1)/rmet(1,1) -rprimd(:,2)*rmet(3,2)/rmet(2,2)
155  width(3)=sqrt(dot_product(vperp,vperp))
156 
157 !Compute the number of cells that will make up the supercell
158  ABI_ALLOCATE(ncells,(natom))
159  mcells=0
160  do iatom=1,natom
161    itypat=typat(iatom)
162    maxrad=radii(npoint(itypat),itypat)
163 !  Compute the lower and higher indices of the supercell
164 !  for this atom
165    do ii=1,3
166      lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
167      highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
168 !    Next coding, still incorrect
169 !    lowest(ii)=floor(xred(ii,iatom)-maxrad/width(ii))-1
170 !    highest(ii)=ceiling(xred(ii,iatom)+maxrad/width(ii))+1
171 !    Old coding incorrect
172 !    lowest(ii)=ceiling(-xred(ii,iatom)-maxrad/width(ii))
173 !    highest(ii)=floor(-xred(ii,iatom)+maxrad/width(ii)+1)
174    end do
175    ncells(iatom)=(highest(1)-lowest(1)+1)* &
176 &   (highest(2)-lowest(2)+1)* &
177 &   (highest(3)-lowest(3)+1)
178 !  DEBUG
179 !  write(std_out,*)' maxrad=',maxrad
180 !  write(std_out,*)' lowest(:)=',lowest(:)
181 !  write(std_out,*)' highest(:)=',highest(:)
182 !  write(std_out,*)' ncells(iatom)=',ncells(iatom)
183 !  ENDDEBUG
184  end do
185  mcells=maxval(ncells(:))
186 
187 !Compute, for each atom, the set of image atoms in the whole supercell
188  ABI_ALLOCATE(xcartcells,(3,mcells,natom))
189  do iatom=1,natom
190    itypat=typat(iatom)
191    maxrad=radii(npoint(itypat),itypat)
192 !  Compute the lower and higher indices of the supercell
193 !  for this atom
194 
195    do ii=1,3
196      lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
197      highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
198    end do
199    icell=0
200    do i1=lowest(1),highest(1)
201      do i2=lowest(2),highest(2)
202        do i3=lowest(3),highest(3)
203          icell=icell+1
204          xcartcells(:,icell,iatom)=xcart(:,iatom)+i1*rprimd(:,1)+i2*rprimd(:,2)+i3*rprimd(:,3)
205        end do
206      end do
207    end do
208  end do
209 
210 !Compute, for each atom, the all-electron pro-atom density
211 !at each point in the primitive cell
212  ABI_ALLOCATE(local_den,(nrx,nry,nrz,natom))
213  ABI_ALLOCATE(step,(2,mpoint))
214  ABI_ALLOCATE(work,(mpoint))
215  ABI_ALLOCATE(yder2,(mpoint))
216  ABI_ALLOCATE(coord1,(3,nrx))
217  coeff1=one/nrx
218  coeff2=one/nry
219  coeff3=one/nrz
220 
221  do iatom=1,natom
222    itypat=typat(iatom)
223    npt=npoint(itypat)
224    maxrad=radii(npt,itypat)
225 !   write(std_out,*)
226 !   write(std_out,'(a,i4)' )' hirsh : accumulating density for atom ',iatom
227 !  write(std_out,*)' ncells(iatom)=',ncells(iatom)
228    do istep=1,npt-1
229      step(1,istep)=radii(istep+1,itypat) - radii(istep,itypat)
230      step(2,istep)=one/step(1,istep)
231    end do
232 !  Approximate first derivative for small radii
233    yp1=(aeden(2,itypat)-aeden(1,itypat))/(radii(2,itypat)-radii(1,itypat))
234    ypn=zero
235    call spline(radii(1:npt,itypat),aeden(1:npt,itypat),npt,yp1,ypn,yder2)
236 
237    local_den(:,:,:,iatom)=zero
238 
239 !  Big loop on the cells
240    do icell=1,ncells(iatom)
241 !    write(std_out,*)' icell=',icell
242      coordat(:)=xcartcells(:,icell,iatom)
243 
244 !    Big loop on the grid points
245      do k1 = 1,nrx
246        coord1(:,k1)=rprimd(:,1)*(k1-1)*coeff1
247      end do
248      do k3 = 1, nrz
249        do k2 = 1, nry
250          coord23_1=rprimd(1,2)*(k2-1)*coeff2+rprimd(1,3)*(k3-1)*coeff3-coordat(1)
251          coord23_2=rprimd(2,2)*(k2-1)*coeff2+rprimd(2,3)*(k3-1)*coeff3-coordat(2)
252          coord23_3=rprimd(3,2)*(k2-1)*coeff2+rprimd(3,3)*(k3-1)*coeff3-coordat(3)
253          do k1 = 1, nrx
254            diff1=coord1(1,k1)+coord23_1
255            diff2=coord1(2,k1)+coord23_2
256            diff3=coord1(3,k1)+coord23_3
257            rr2=diff1**2+diff2**2+diff3**2
258            if(rr2<maxrad**2)then
259 
260              rr=sqrt(rr2)
261 !            Find the index of the radius by bissection
262              if (rr < radii(1,itypat)) then
263 !              Linear extrapolation
264                den=aeden(1,itypat)+(rr-radii(1,itypat))/(radii(2,itypat)-radii(1,itypat))&
265 &               *(aeden(2,itypat)-aeden(1,itypat))
266              else
267 !              Use the spline interpolation
268 !              Find the index of the radius by bissection
269                inmin=1
270                inmax=npt
271                igrid=1
272                do
273                  if(inmax-inmin==1)exit
274                  igrid=(inmin+inmax)/2
275                  if(rr>=radii(igrid,itypat))then
276                    inmin=igrid
277                  else
278                    inmax=igrid
279                  end if
280                end do
281                igrid=inmin
282 !              write(std_out,*)' igrid',igrid
283 
284                hh=step(1,igrid)
285                h_inv=step(2,igrid)
286                aa= (radii(igrid+1,itypat)-rr)*h_inv
287                bb= (rr-radii(igrid,itypat))*h_inv
288                den = aa*aeden(igrid,itypat) + bb*aeden(igrid+1,itypat)  &
289 &               +( (aa*aa*aa-aa)*yder2(igrid)         &
290 &               +(bb*bb*bb-bb)*yder2(igrid+1) ) *hh*hh*sixth
291              end if ! Select small radius or spline
292 
293              local_den(k1,k2,k3,iatom)=local_den(k1,k2,k3,iatom)+den
294            end if ! dist2<maxrad
295 
296          end do ! k1
297        end do ! k2
298      end do ! k3
299 
300    end do ! icell
301  end do ! iatom
302 
303 !Compute, the total all-electron density at each point in the primitive cell
304  ABI_ALLOCATE(sum_den,(nrx,nry,nrz))
305  sum_den(:,:,:)=zero
306  do iatom=1,natom
307    sum_den(:,:,:)=sum_den(:,:,:)+local_den(:,:,:,iatom)
308  end do
309 
310 !DEBUG
311 !do k3=1,nrz
312 !do k2=1,nry
313 !do k1=1,nrx
314 !write(std_out,'(3i4,3es16.6)' )k1,k2,k3,local_den(k1,k2,k3,1:2),sum_den(k1,k2,k3)
315 !end do
316 !end do
317 !end do
318 !write(std_out,*)' hirsh : before accumulate the integral of the density'
319 !ENDDEBUG
320 
321 
322 !Accumulate the integral of the density, to get Hirshfeld charges
323 !There is a minus sign because the electron has a negative charge
324  ngoodpoints = 0
325  hcharge(:)=zero
326  hweight(:)=zero
327  do k3=1,nrz
328    do k2=1,nry
329      do k1=1,nrx
330 !      Use minimal_den in order to avoid divide by zero
331        if (abs(sum_den(k1,k2,k3)) > minimal_den) then
332          ngoodpoints = ngoodpoints+1
333          factor=grid_den(k1,k2,k3)/(sum_den(k1,k2,k3)+minimal_den)
334          do iatom=1,natom
335            hden(iatom)=hden(iatom)+local_den(k1,k2,k3,iatom)
336            hcharge(iatom)=hcharge(iatom)-local_den(k1,k2,k3,iatom)*factor
337            hweight(iatom)=hweight(iatom)+local_den(k1,k2,k3,iatom)/(sum_den(k1,k2,k3)+minimal_den)
338          end do
339        end if
340      end do
341    end do
342  end do
343 
344 !DEBUG
345 !do iatom=1,natom
346 !write(std_out,'(i9,3es17.6)' )iatom,hden(iatom),hcharge(iatom),hweight(iatom)
347 !end do
348 !ENDDEBUG
349 
350  hcharge(:)=hcharge(:)*ucvol/dble(nfftot)
351  hweight(:)=hweight(:)/dble(nfftot)
352 
353 !Check on the total charge
354  total_zion=sum(zion(typat(1:natom)))
355  total_charge=sum(hcharge(1:natom))
356  total_weight=sum(hweight(1:natom))
357 
358 !DEBUG
359 !write(std_out,*)' ngoodpoints = ', ngoodpoints, ' out of ', nfftot
360 !write(std_out,*)' total_weight=',total_weight
361 !write(std_out,*)' total_weight=',total_weight
362 !ENDDEBUG
363 
364 !Output
365  if (prtcharge == 1) then
366      write(std_out,*)
367      write(std_out,*)'    Hirshfeld analysis'
368      write(std_out,*)'    Atom       Zion       Electron  Charge       Net charge '
369      write(std_out,*)
370      do iatom=1,natom
371        write(std_out,'(i9,3es17.6)' )&
372 &       iatom,zion(typat(iatom)),hcharge(iatom),hcharge(iatom)+zion(typat(iatom))
373      end do
374      write(std_out,*)
375      write(std_out,'(a,3es17.6)')'    Total',total_zion,total_charge,total_charge+total_zion
376      write(std_out,*)
377  end if
378 
379  ABI_DEALLOCATE(coord1)
380  ABI_DEALLOCATE(local_den)
381  ABI_DEALLOCATE(ncells)
382  ABI_DEALLOCATE(step)
383  ABI_DEALLOCATE(sum_den)
384  ABI_DEALLOCATE(work)
385  ABI_DEALLOCATE(xcartcells)
386  ABI_DEALLOCATE(xred)
387  ABI_DEALLOCATE(yder2)
388 
389 end subroutine dens_hirsh