TABLE OF CONTENTS
ABINIT/atomden [ Functions ]
NAME
atomden
FUNCTION
Construct atomic proto-bulk density (i.e. the superposed density
from neutral, isolated atoms at the bulk atomic positions).
This is useful if one wants to construct the bonding density:
rho^{bnd} = rho^{bulk}(r)
- \sum_{\alpha}\rho^{atm}_{\alpha}(r-R_{\alpha})
Where rho^{bulk} is the bulk density, rho^{atm} the atomic density
and the index \alpha sums over all atoms. the R_{\alpha} are the
atomic positions in the bulk. This routine calculates the sum over
rho^{atm}_{\alpha}(r-R_{\alpha}) on a grid.
Units are atomic.
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (SM,VR,FJ,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 .
INPUTS
calctype : type of calculation
'replace' zero the input/output density array
'add' add to the input/output density array
natom : number of atoms in cell
ntypat : number of different types of atoms in cell
typat(natom) : type of each atom
ngrid : number of gridpoints
r_vec_grid(3,ngrid) : real (non-reduced) coordinates for grid points
rho(ngrid) : input/output density array
a(3),b(3),c(3) : real-space basis vectors
atom_pos(3,natom) : reduced coordinates for atomic positions
natomgr(ntypat) : number of gridpoints for each atomic density grid
natomgrmax : max(natomgr(ntypat))
atomrgrid(natomgrmax,ntypat)
density(natomgrmax,ntypat)
OUTPUT
rho(ngrid) : input/output density array
SIDE EFFECTS
NOTES
There are two ways to compile the proto density in real space for a solid. One alternative is that the density is calculated for an extended grid encompassing the sphere of points around one atom, and the results are folded back into the unit cell. On the other hand one can, around each grid point, identify the number of atoms in a sphere equivalent to the length of the radial grid for each type of atom. The second approach, with some modification, is taken here. The numer of atoms in a supercell cell are listed such that the supercell encompasses the atoms which could contribute to any point in the grid. That list is kept and cycled through, to avoid recalculating it at each point. Note that the density calculated from the atom is the spherical average, since there is no preferred direction without any external field (and it's simpler)
PARENTS
read_atomden
CHILDREN
sort_dp,spline,splint,wrtout,xmpi_barrier,xmpi_sum_master
SOURCE
75 #if defined HAVE_CONFIG_H 76 #include "config.h" 77 #endif 78 79 #include "abi_common.h" 80 81 subroutine atomden(MPI_enreg,natom,ntypat,typat,ngrid,r_vec_grid,rho,a,b,c,atom_pos, & 82 & natomgr,natomgrmax,atomrgrid,density,prtvol,calctype) 83 84 use defs_basis 85 use defs_abitypes 86 use m_profiling_abi 87 use m_xmpi 88 use m_splines 89 use m_sort 90 91 !This section has been created automatically by the script Abilint (TD). 92 !Do not modify the following lines by hand. 93 #undef ABI_FUNC 94 #define ABI_FUNC 'atomden' 95 use interfaces_14_hidewrite 96 !End of the abilint section 97 98 implicit none 99 100 !Arguments ------------------------------------ 101 !scalars 102 integer,intent(in) :: natom,ntypat,ngrid,natomgrmax,prtvol 103 character(len=7),intent(in) :: calctype 104 !arrays 105 type(MPI_type),intent(in) :: MPI_enreg 106 integer,intent(in) :: typat(natom),natomgr(ntypat) 107 real(dp),intent(in) :: r_vec_grid(3,ngrid),a(3),b(3),c(3) 108 real(dp),intent(in) :: atom_pos(3,natom),atomrgrid(natomgrmax,ntypat) 109 real(dp),intent(in) :: density(natomgrmax,ntypat) 110 real(dp),intent(inout) :: rho(ngrid) 111 112 !Local variables------------------------------- 113 !scalars 114 character(len=500) :: message 115 integer :: cnt,delta,i,l,m,n,iatom,itypat,igrid,n_cells,n_grid_p 116 integer :: ierr,spaceComm,nprocs,master,rank,remainder 117 real(dp) :: a_norm,b_norm,c_norm 118 real(dp) :: r_max,R_sphere_max,dp_dummy,ybcbeg,ybcend 119 !arrays 120 integer :: n_equiv_atoms(ntypat),grid_index(ngrid) 121 integer :: my_start_equiv_atoms(ntypat) 122 integer :: my_end_equiv_atoms(ntypat) 123 integer :: l_min(ntypat),m_min(ntypat),n_min(ntypat) 124 integer :: l_max(ntypat),m_max(ntypat),n_max(ntypat) 125 real(dp) :: center(3),dp_vec_dummy(3),delta_a(3),delta_b(3),delta_c(3) 126 real(dp) :: r_atom(3),grid_distances(ngrid) 127 integer, allocatable :: new_index(:),i_1d_dummy(:) 128 real(dp),allocatable :: equiv_atom_dist(:,:),equiv_atom_pos(:,:,:),rho_temp(:,:) 129 real(dp),allocatable :: dp_1d_dummy(:),dp_2d_dummy(:,:),ypp(:) 130 real(dp),allocatable :: x_fit(:),y_fit(:) 131 132 133 ! ************************************************************************ 134 135 !initialise and check parallel execution 136 spaceComm=MPI_enreg%comm_cell 137 nprocs=xmpi_comm_size(spaceComm) 138 rank=MPI_enreg%me_kpt 139 140 master=0 141 142 !initialise variables and vectors 143 a_norm = sqrt(dot_product(a,a)) 144 b_norm = sqrt(dot_product(b,b)) 145 c_norm = sqrt(dot_product(c,c)) 146 center = (a+b+c)*half 147 dp_dummy = dot_product(a,b)/(b_norm*b_norm) 148 dp_vec_dummy = dp_dummy*b 149 delta_a = a - dp_vec_dummy 150 dp_dummy = dot_product(b,a)/(a_norm*a_norm) 151 dp_vec_dummy = dp_dummy*a 152 delta_b = b - dp_vec_dummy 153 dp_dummy = dot_product(c,(a+b))/(dot_product((a+b),(a+b))) 154 dp_vec_dummy = dp_dummy*(a+b) 155 delta_c = c - dp_vec_dummy 156 ABI_ALLOCATE(rho_temp,(ngrid,ntypat)) 157 rho_temp = zero 158 159 !write(std_out,*) '*** --- In atomden --- ***' 160 !write(std_out,*) ' a_norm:',a_norm,' b_norm:',b_norm,' c_norm:',c_norm 161 !write(std_out,*) 'delta_a:',delta_a,'delta_b:',delta_b,'delta_c:',delta_c 162 !write(std_out,*) ' center:',center 163 164 !Find supercell which will contain all possible contributions 165 !for all atoms, and enumerate positions for all atoms 166 !TODO list of atoms can be "pruned", i.e identify all atoms 167 !that can't possibly contribute and remove from list. 168 !Should be most important for very oblique cells 169 do itypat=1,ntypat 170 R_sphere_max = atomrgrid(natomgr(itypat),itypat) 171 l_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_a,delta_a))) 172 l_max(itypat) = -l_min(itypat) 173 m_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_b,delta_b))) 174 m_max(itypat) = -m_min(itypat) 175 n_min(itypat) = -ceiling(R_sphere_max/sqrt(dot_product(delta_c,delta_c))) 176 n_max(itypat) = -n_min(itypat) 177 n_cells = (l_max(itypat)-l_min(itypat)+1) & 178 & *(m_max(itypat)-m_min(itypat)+1) & 179 & *(n_max(itypat)-n_min(itypat)+1) 180 n_equiv_atoms(itypat) = 0 181 do iatom=1,natom 182 if (typat(iatom)==itypat) then 183 n_equiv_atoms(itypat) = n_equiv_atoms(itypat) + n_cells 184 end if ! if type=itypat 185 end do ! number of atoms per cell 186 if ((rank==master).and.(prtvol>9)) then 187 write(message,'(a)') '*** --- In atomden --- find box ***' 188 call wrtout(std_out,message,'COLL') 189 write(message,'(a,I4)') ' itypat:',itypat 190 call wrtout(std_out,message,'COLL') 191 write(message,'(2(a,I4))') ' l_min:',l_min(itypat),' l_max:',l_max(itypat) 192 call wrtout(std_out,message,'COLL') 193 write(message,'(2(a,I4))') ' m_min:',m_min(itypat),' m_max:',m_max(itypat) 194 call wrtout(std_out,message,'COLL') 195 write(message,'(2(a,I4))') ' n_min:',n_min(itypat),' n_max:',n_max(itypat) 196 call wrtout(std_out,message,'COLL') 197 write(message,'(2(a,I4))') ' n_equiv_atoms:',n_equiv_atoms(itypat) 198 call wrtout(std_out,message,'COLL') 199 end if 200 end do !atom type 201 202 !allocate arrays 203 n = maxval(n_equiv_atoms) 204 ABI_ALLOCATE(equiv_atom_pos,(3,n,ntypat)) 205 ABI_ALLOCATE(equiv_atom_dist,(n,ntypat)) 206 equiv_atom_pos = zero 207 equiv_atom_dist = zero 208 209 !Find positions and distance of atoms from center of cell 210 do itypat=1,ntypat 211 i = 1 212 do l=l_min(itypat),l_max(itypat) 213 do m=m_min(itypat),m_max(itypat) 214 do n=n_min(itypat),n_max(itypat) 215 do iatom=1,natom 216 if (typat(iatom)==itypat) then 217 if (i>n_equiv_atoms(itypat)) then 218 MSG_ERROR('atomden: i>n_equiv_atoms') 219 end if 220 equiv_atom_pos(:,i,itypat) = (atom_pos(1,iatom)+dble(l))*a & 221 & + (atom_pos(2,iatom)+dble(m))*b & 222 & + (atom_pos(3,iatom)+dble(n))*c 223 dp_vec_dummy = equiv_atom_pos(:,i,itypat)-center 224 equiv_atom_dist(i,itypat) = & 225 & sqrt(dot_product(dp_vec_dummy,dp_vec_dummy)) 226 i = i + 1 227 end if 228 end do 229 end do !n 230 end do !m 231 end do !l 232 ! write(std_out,*) '*** --- In atomden --- find equiv ***' 233 ! write(std_out,*) ' itypat:',itypat 234 ! write(std_out,*) ' equiv_atom_pos:' 235 ! write(std_out,*) equiv_atom_pos(:,:,itypat) 236 ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat) 237 end do !atom type 238 239 !Sort the atoms after distance so that the density from the ones 240 !furthest away can be added first. This is to prevent truncation error. 241 do itypat=1,ntypat 242 n = n_equiv_atoms(itypat) 243 ABI_ALLOCATE(dp_1d_dummy,(n)) 244 ABI_ALLOCATE(new_index,(n)) 245 ABI_ALLOCATE(dp_2d_dummy,(3,n)) 246 dp_1d_dummy = equiv_atom_dist(1:n,itypat) 247 dp_2d_dummy = equiv_atom_pos(1:3,1:n,itypat) 248 do i=1,n 249 new_index(i) = i 250 end do 251 call sort_dp(n,dp_1d_dummy,new_index,tol14) 252 do i=1,n 253 ! write(std_out,*) i,' -> ',new_index(i) 254 equiv_atom_pos(1:3,n+1-i,itypat) = dp_2d_dummy(1:3,new_index(i)) 255 equiv_atom_dist(1:n,itypat) = dp_1d_dummy 256 end do 257 ABI_DEALLOCATE(dp_1d_dummy) 258 ABI_DEALLOCATE(new_index) 259 ABI_DEALLOCATE(dp_2d_dummy) 260 ! write(std_out,*) '*** --- In atomden --- sorting atoms ***' 261 ! write(std_out,*) ' itypat:',itypat 262 ! write(std_out,*) ' equiv_atom_pos:' 263 ! write(std_out,*) equiv_atom_pos(:,:,itypat) 264 ! write(std_out,*) ' equiv_atom_dist:',equiv_atom_dist(:,itypat) 265 end do ! atom type 266 267 !Divide the work in case of parallel execution 268 if (nprocs==1) then ! Make sure everything runs with one proc 269 if (prtvol>9) then 270 write(message,'(a)') ' In atomden - number of processors: 1' 271 call wrtout(std_out,message,'COLL') 272 write(message,'(a)') ' Calculation of proto-atomic density done in serial' 273 call wrtout(std_out,message,'COLL') 274 end if 275 do itypat=1,ntypat 276 if (prtvol>9) then 277 write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat) 278 call wrtout(std_out,message,'COLL') 279 end if 280 my_start_equiv_atoms(itypat) = 1 281 my_end_equiv_atoms(itypat) = n_equiv_atoms(itypat) 282 end do 283 else 284 if (rank==master.and.prtvol>9) then 285 write(message,'(a,I5)') ' In atomden - number of processors:',nprocs 286 call wrtout(std_out,message,'COLL') 287 write(message,'(a)') ' Calculation of proto-atomic density done in parallel' 288 call wrtout(std_out,message,'COLL') 289 end if 290 do itypat=1,ntypat 291 if (rank==master.and.prtvol>9) then 292 write(message,'(a,I6)') ' Number of equivalent atoms:',n_equiv_atoms(itypat) 293 call wrtout(std_out,message,'COLL') 294 end if 295 ! Divide the atoms among the processors by shuffling indices 296 delta = int(floor(real(n_equiv_atoms(itypat))/real(nprocs))) 297 remainder = n_equiv_atoms(itypat)-nprocs*delta 298 my_start_equiv_atoms(itypat) = 1+rank*delta 299 my_end_equiv_atoms(itypat) = (rank+1)*delta 300 ! Divide the remainder points among the processors 301 ! by shuffling indices 302 if ((rank+1)>remainder) then 303 my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + remainder 304 my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + remainder 305 else 306 my_start_equiv_atoms(itypat) = my_start_equiv_atoms(itypat) + rank 307 my_end_equiv_atoms(itypat) = my_end_equiv_atoms(itypat) + rank + 1 308 end if 309 if (prtvol>9) then 310 write(message,'(a,I3)') ' For atom type: ',itypat 311 call wrtout(std_out,message,'PERS') 312 ! write(message,'(a,I6)') ' I''ll take atoms from: ',my_start_equiv_atoms(itypat) 313 ! call wrtout(std_out,message,'PERS') 314 ! write(message,'(a,I6)') ' total for me: ',my_end_equiv_atoms(itypat) 315 ! call wrtout(std_out,message,'PERS') 316 write(message,'(a,I6)') ' total for me: ', & 317 & my_end_equiv_atoms(itypat)+1-my_start_equiv_atoms(itypat) 318 call wrtout(std_out,message,'PERS') 319 end if 320 end do 321 end if 322 323 !Loop over types of atoms and equivalent atoms and 324 !interpolate density onto grid 325 do itypat=1,ntypat 326 ! do iatom=my_start_equiv_atoms(itypat),my_end_equiv_atoms(itypat) 327 328 cnt = 0 329 iatom = rank+1 - nprocs 330 ! Parallel execution of loop 331 do 332 cnt = cnt + 1 333 iatom = iatom + nprocs 334 if (iatom>n_equiv_atoms(itypat)) exit ! Exit if index is too large 335 336 if (mod(cnt,100)==0.and.prtvol>0) then 337 write(message,'(2(a,I6))') ' atoms so far',cnt,' of: ',n_equiv_atoms(itypat)/nprocs 338 call wrtout(std_out,message,'PERS') 339 end if 340 341 r_max = atomrgrid(natomgr(itypat),itypat) 342 r_atom = equiv_atom_pos(:,iatom,itypat) 343 344 ! Set up an array with the gridpoint distances 345 i = 1 346 grid_distances = zero 347 grid_index = 0 348 do igrid=1,ngrid 349 dp_vec_dummy(:) = r_vec_grid(:,igrid) - r_atom(:) 350 dp_dummy = sqrt(dot_product(dp_vec_dummy,dp_vec_dummy)) 351 if (dp_dummy <= r_max) then 352 grid_distances(i) = dp_dummy 353 grid_index(i) = igrid 354 i = i + 1 355 else 356 cycle ! cycle if point is too far away 357 end if 358 end do 359 n_grid_p = i - 1 360 361 if (n_grid_p==0) cycle ! Cycle if no point needs 362 ! to be interpolated 363 364 ! Sort points to be interpolated in ascending order 365 ABI_ALLOCATE(dp_1d_dummy,(n_grid_p)) 366 ABI_ALLOCATE(new_index,(n_grid_p)) 367 ABI_ALLOCATE(i_1d_dummy,(n_grid_p)) 368 dp_1d_dummy = grid_distances(1:n_grid_p) 369 do i=1,n_grid_p 370 new_index(i) = i 371 end do 372 call sort_dp(n_grid_p,dp_1d_dummy,new_index,tol16) 373 grid_distances(1:n_grid_p) = dp_1d_dummy 374 i_1d_dummy = grid_index(1:n_grid_p) 375 do i=1,n_grid_p 376 ! write(std_out,*) i_1d_dummy(i),' -> ',i_1d_dummy(new_index(i)) 377 grid_index(i) = i_1d_dummy(new_index(i)) 378 end do 379 ABI_DEALLOCATE(dp_1d_dummy) 380 ABI_DEALLOCATE(new_index) 381 ABI_DEALLOCATE(i_1d_dummy) 382 383 ! Interpolate density onto all grid points 384 ABI_ALLOCATE(ypp,(natomgr(itypat))) 385 ABI_ALLOCATE(x_fit,(n_grid_p)) 386 ABI_ALLOCATE(y_fit,(n_grid_p)) 387 ypp = zero; y_fit = zero 388 ybcbeg = zero; ybcend = zero 389 x_fit = grid_distances(1:n_grid_p) 390 call spline(atomrgrid(1:natomgr(itypat),itypat), & 391 & density(1:natomgr(itypat),itypat), & 392 & natomgr(itypat),ybcbeg,ybcend,ypp) 393 call splint(natomgr(itypat),atomrgrid(1:natomgr(itypat),itypat), & 394 & density(1:natomgr(itypat),itypat),ypp,n_grid_p, & 395 & x_fit,y_fit) 396 397 ! Save the interpolated points to grid 398 do i=1,n_grid_p 399 rho_temp(grid_index(i),itypat) = rho_temp(grid_index(i),itypat) + y_fit(i) 400 end do 401 ABI_DEALLOCATE(ypp) 402 ABI_DEALLOCATE(x_fit) 403 ABI_DEALLOCATE(y_fit) 404 405 end do ! n equiv atoms 406 end do ! type of atom 407 408 !Collect all contributions to rho_temp if 409 !we are running in parallel 410 if (nprocs>1) then 411 call xmpi_barrier(spaceComm) 412 call xmpi_sum_master(rho_temp,master,spaceComm,ierr) 413 call xmpi_barrier(spaceComm) 414 if (prtvol>9) then 415 write(message,'(a)') ' In atomden - contributions to rho_temp collected' 416 call wrtout(std_out,message,'PERS') 417 end if 418 end if 419 420 !Now rho_temp contains the atomic protodensity for each atom. 421 !Check whether this is to replace or be added to the input/output array 422 !and sum up contributions 423 if (trim(calctype)=='replace') rho = zero 424 do itypat=1,ntypat 425 rho(:) = rho(:) + rho_temp(:,itypat) 426 end do 427 428 !deallocations 429 if (allocated(rho_temp)) then 430 ABI_DEALLOCATE(rho_temp) 431 end if 432 if (allocated(equiv_atom_pos)) then 433 ABI_DEALLOCATE(equiv_atom_pos) 434 end if 435 if (allocated(equiv_atom_dist)) then 436 ABI_DEALLOCATE(equiv_atom_dist) 437 end if 438 !if (allocated()) deallocate() 439 440 return 441 442 end subroutine atomden