TABLE OF CONTENTS
ABINIT/m_dens [ 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