1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 module atomlist 9 10 use precision 11 use alloc 12 use parallel, only: IOnode 13 use atmfuncs, only: nofis, nkbfis, izofis, massfis, 14 $ rcut, atmpopfio, zvalfis, floating 15 use atm_types, only: species 16 use siesta_geom, only: na_u, na_s, xa, isa, xa_last 17 implicit none 18 19 private 20 public :: initatomlists, superc 21 public :: superx ! for backwards compatibility 22 public :: reset_atomlists 23 24! 25! Instead of "generic" na, no, and nokb, we use: 26! 27! For "supercell" (intended for k-point calcs) 28 integer, save, public :: no_s ! Number of orbitals 29 integer, save, public :: nokb_s ! Number of KB projs 30 31! Same for "unit", or "real" cell: 32 integer, save, public :: no_u, nokb_u 33 integer, save, public :: no_l=1 ! Local to node 34 35! Here 'na' is a generic number. It could be na_u or na_s, depending 36! on whether we need a supercell or not. 37 38C character*2 elem(na) : Element name of each atom. 39C integer lasto(0:na) : Position of last orbital of each atom 40C integer lastkb(0:na) : Position of last KB proj. of each atom 41C integer iza(na) : Atomic number of each atom 42C real*8 amass(na) : Atomic mass of each atom 43C real*8 qa(na) : Neutral atom charge of each atom 44 45 character(len=2), pointer, save, public :: elem(:) => null() 46! elem will contain element names, so is 2 chars in length 47 48 integer, pointer, save, public :: iza(:) => null() 49 integer, pointer, save, public :: lasto(:) => null() 50 integer, pointer, save, public :: lastkb(:) => null() 51 real(dp), pointer, save, public :: amass(:) => null() 52 real(dp), pointer, save, public :: qa(:) => null() 53 integer, pointer, save, public :: indxua(:) => null() 54 55 ! This array depends on the actual geometry, so 56 ! it does not properly belong here. It is used only 57 ! for communication between nlefsm and hsparse. 58 ! For safety, it is initialized in each invokation 59 ! of hsparse. 60 logical, pointer, save, public :: in_kb_orb_u_range(:) => null() 61 62! Index of equivalent atom in "u" cell 63 real(dp), save, public :: rmaxv ! Max cutoff for Vna 64 real(dp), save, public :: rmaxo ! Max cuoff for at. orb. 65 real(dp), save, public :: rmaxkb ! Max cuoff for KB proj. 66 real(dp), save, public :: rmaxdftu ! Max cuoff for DFTU proj. 67 68 real(dp), save, public :: qtot ! Total number of elect. 69 real(dp), save, public :: qtots(2) ! Total number of electrons per spin 70 71 real(dp), save, public :: zvaltot 72 ! Total number of pseudoprotons 73 ! (excluding those of ghost atoms) 74 75 76 integer, pointer, save, public :: iaorb(:) => null() 77 ! Atomic index of each orbital 78 integer, pointer, save, public :: iphorb(:) => null() 79 ! Orbital index of each orbital in its atom 80 real(dp), pointer, save, public :: Datm(:) => null() 81 ! Neutral atom charge 82 ! of each orbital 83 real(dp), pointer, save, public :: rco(:) => null() 84 ! Cutoff radius of each orbital 85 86 integer, pointer, save, public :: indxuo(:) => null() 87 ! Index of equivalent orbital in "u" cell 88 89 integer, pointer, save, public :: iakb(:) => null() 90! Atomic index of each KB projector 91 integer, pointer, save, public :: iphKB(:) => null() 92! Index of each KB projector in its atom (negative) 93 real(dp), pointer, save, public :: rckb(:) => null() 94! Cutoff radius of each KB projector 95! 96 97 CONTAINS 98 99!======================================================= 100 subroutine initatomlists() 101 102 use atm_types, only: species_info 103 use radial, only: rad_func 104 use dftu_specs, only: switch_dftu 105 106C Routine to initialize the atomic lists. 107C 108 integer ia, io, is, nkba, noa, nol, nokbl, ioa, ikb 109 type(species_info), pointer :: spp 110 type(rad_func), pointer :: pp 111 112 call re_alloc( indxua, 1, na_u, 'indxua', 'atomlist' ) 113 call re_alloc( lastkb, 0, na_u, 'lastkb', 'atomlist' ) 114 call re_alloc( lasto, 0, na_u, 'lasto', 'atomlist' ) 115 call re_alloc( qa, 1, na_u, 'qa', 'atomlist' ) 116 call re_alloc(xa_last,1,3,1,na_u,'xa_last','atomlist') 117 call re_alloc( amass, 1, na_u, 'amass', 'atomlist' ) 118 call re_alloc( in_kb_orb_u_range, 1, na_u, 'in_kb', 'atomlist' ) 119! 120! Find number of orbitals and KB projectors in cell 121! 122 no_u = 0 123 nokb_u = 0 124 do ia = 1,na_u 125 is = isa(ia) 126 noa = species(is)%norbs 127 nkba = species(is)%nprojs 128 no_u = no_u + noa 129 nokb_u = nokb_u + nkba 130 enddo 131 132 na_s = na_u 133 no_s = no_u 134 nokb_s = nokb_u 135 136 call re_alloc( iaorb, 1, no_u, 'iaorb', 'atomlist' ) 137 call re_alloc( indxuo, 1, no_u, 'indxuo', 'atomlist' ) 138 call re_alloc( iphorb, 1, no_u, 'iphorb', 'atomlist' ) 139 call re_alloc( Datm, 1, no_u, 'Datm', 'atomlist' ) 140 call re_alloc( rco, 1, no_u, 'rco', 'atomlist' ) 141 call re_alloc( iaKB, 1, nokb_u, 'iaKB', 'atomlist' ) 142 call re_alloc( iphKB, 1, nokb_u, 'iphKB', 'atomlist' ) 143 call re_alloc( rckb, 1, nokb_u, 'rckb', 'atomlist' ) 144 145c Initialize atomic lists 146 nol = 0 147 nokbl = 0 148 qtot = 0._dp 149 rmaxv = 0._dp 150 rmaxo = 0._dp 151 rmaxkb = 0._dp 152 rmaxdftu = 0._dp 153 lasto(0) = 0 154 lastkb(0) = 0 155 zvaltot = 0.0_dp 156 do ia = 1,na_u 157 is = isa(ia) 158 if (.not. floating(is)) then 159 zvaltot = zvaltot + zvalfis(is) 160 endif 161 noa = nofis(is) 162 nkba = nkbfis(is) 163 lasto(ia) = lasto(ia-1) + noa 164 lastkb(ia) = lastkb(ia-1) + nkba 165 rmaxv = max( rmaxv, rcut(is,0) ) 166 iza(ia) = izofis(is) 167 amass(ia) = massfis(is) 168 qa(ia) = 0.0_dp 169 do io = 1,noa 170 nol = nol + 1 171 rmaxo = max( rmaxo, rcut(is,io) ) 172 iaorb(nol) = ia 173 iphorb(nol) = io 174 Datm(nol) = atmpopfio(is,io) 175 qa(ia) = qa(ia) + Datm(nol) 176 qtot = qtot + Datm(nol) 177 enddo 178 do io = 1,nkba 179 nokbl = nokbl + 1 180 rmaxkb = max( rmaxkb, rcut(is,-io) ) 181 iaKB(nokbl) = ia 182 iphKB(nokbl) = -io 183 enddo 184 if( switch_dftu ) then 185 spp => species(is) 186 do io = 1, spp%n_pjdftunl 187 pp => spp%pjdftu(io) 188 rmaxdftu = max( rmaxdftu, pp%cutoff ) 189 enddo 190 endif 191 enddo 192 193! Find rco and rckb ............................. 194 195 do ia = 1,na_u 196 is = isa(ia) 197 do io = lasto(ia-1)+1,lasto(ia) 198 ioa = iphorb(io) 199 rco(io) = rcut(is,ioa) 200 enddo 201 do ikb = lastkb(ia-1)+1,lastkb(ia) 202 ioa = iphKB(ikb) 203 rckb(ikb) = rcut(is,ioa) 204 enddo 205 enddo 206 207 if (IOnode) 208 $ write(6,'(/a,3(1x,i5))') 209 $ 'initatomlists: Number of atoms, orbitals, and projectors: ', 210 $ na_u, no_u, nokb_u 211 212 end subroutine initatomlists 213 214 subroutine reset_atomlists() 215 216 use alloc, only: de_alloc 217 218 call de_alloc( indxua, 'indxua', 'atomlist' ) 219 call de_alloc( isa, 'isa', 'atomlist') 220 call de_alloc( iza, 'iza', 'atomlist') 221 call de_alloc( lastkb, 'lastkb', 'atomlist' ) 222 call de_alloc( lasto, 'lasto', 'atomlist' ) 223 call de_alloc( qa, 'qa', 'atomlist' ) 224 call de_alloc(xa_last,'xa_last','atomlist') 225 call de_alloc( amass, 'amass', 'atomlist' ) 226 call de_alloc( in_kb_orb_u_range, 'in_kb', 'atomlist' ) 227 228 call de_alloc( iaorb, 'iaorb', 'superc' ) 229 call de_alloc( indxuo, 'indxuo', 'superc' ) 230 call de_alloc( iphorb, 'iphorb', 'superc' ) 231 call de_alloc( Datm, 'Datm', 'superc' ) 232 call de_alloc( rco, 'rco', 'superc' ) 233 call de_alloc( iaKB, 'iaKB', 'superc' ) 234 call de_alloc( iphKB, 'iphKB', 'superc' ) 235 call de_alloc( rckb, 'rckb', 'superc' ) 236 237 end subroutine reset_atomlists 238 239 subroutine superc( ucell, scell, nsc) 240 241C Finds the supercell required to avoid multiple image overlaps, 242C and expands arrays from unit cell to supercell 243C Written by J.M.Soler. August 1998. 244! Rewritten Alberto Garcia, May 2000. 245 246 implicit none 247 integer, intent(in) :: nsc(3) ! Diagonal elements of mscell 248 real(dp), intent(in) :: ucell(3,3) ! Unit cell vectors 249 real(dp), intent(out) :: scell(3,3) ! Supercell vectors 250 251C Internal variables 252 integer ia, io, iua, iuo, ja, ncells, 253 $ na, no, nokb 254 255! 256! Find number of cells, atoms and orbitals in supercell 257 ncells = nsc(1) * nsc(2) * nsc(3) 258 na = na_u * ncells 259 no = no_u * ncells 260 nokb = nokb_u * ncells 261 262! 263! Reallocate arrays if needed 264! 265 if (na.gt.na_s) then 266 call re_alloc( indxua, 1, na, 'indxua', 'atomlist', .true. ) 267 call re_alloc( isa, 1, na, 'isa', 'atomlist', .true. ) 268 call re_alloc( iza, 1, na, 'iza', 'atomlist', .true. ) 269 call re_alloc( lastkb, 0, na, 'lastkb', 'atomlist', .true. ) 270 call re_alloc( lasto, 0, na, 'lasto', 'atomlist', .true. ) 271 call re_alloc( in_kb_orb_u_range, 1, na, 272 $ 'in_kb', 'atomlist', .true. ) 273 call re_alloc( qa, 1, na, 'qa', 'atomlist', .true. ) 274 call re_alloc( xa, 1, 3, 1, na, 'xa', 'atomlist', .true. ) 275 call re_alloc(xa_last, 1,3, 1,na, 'xa_last', 'superc', 276 & copy=.true. ) 277 endif 278 279 na_s = na 280 281C Find supercell vectors and atomic coordinates in supercell 282 call superx( ucell, nsc, na_u, na_s, xa, scell ) 283 284C Find indxua and expand isa, iza, lasto and lastkb to supercell 285 do ia = 1,na_s 286 ja = mod(ia-1,na_u) + 1 287 indxua(ia) = ja 288 isa(ia) = isa(ja) 289 iza(ia) = iza(ja) 290 lasto(ia) = lasto(ia-1) + lasto(ja) - lasto(ja-1) 291 lastkb(ia) = lastkb(ia-1) + lastkb(ja) - lastkb(ja-1) 292 enddo 293 294! Reallocate orbital arrays 295 296 if (no.gt.no_s) then 297 call re_alloc(iaorb, 1,no, routine='superc',copy=.true.) 298 call re_alloc(indxuo, 1,no, routine='superc',copy=.true.) 299 call re_alloc(iphorb, 1,no, routine='superc',copy=.true.) 300 call re_alloc(Datm, 1,no, routine='superc',copy=.true.) 301 call re_alloc(rco, 1,no, routine='superc',copy=.true.) 302 endif 303 304 no_s = no 305 306C Find indxuo and expand iaorb, iphorb, and rco 307 do io = 1,no_s 308 indxuo(io) = mod(io-1,no_u) + 1 309 enddo 310 do ia = 1,na_s 311 do io = lasto(ia-1)+1,lasto(ia) 312 iuo = indxuo(io) 313 iaorb(io) = ia 314 iphorb(io) = iphorb(iuo) 315 rco(io) = rco(iuo) 316 enddo 317 enddo 318 319! Reallocate projector arrays 320 321 if (nokb .gt. nokb_s) then 322 call re_alloc(iaKB, 1,nokb, routine='superc',copy=.true.) 323 call re_alloc(iphKB, 1,nokb, routine='superc',copy=.true.) 324 call re_alloc(rckb, 1,nokb, routine='superc',copy=.true.) 325 endif 326 327 nokb_s = nokb 328 329C Expand iakb and iphKB and rckb 330 331 do ia = 1,na_s 332 iua = indxua(ia) 333 iuo = lastkb(iua-1) 334 do io = lastkb(ia-1)+1,lastkb(ia) 335 iuo = iuo + 1 336 iakb(io) = ia 337 iphKB(io) = iphKB(iuo) 338 rckb(io) = rckb(iuo) 339 enddo 340 enddo 341 342 if (IOnode .and. ncells.gt.1) then 343 write(6,'(/,a,i6,a,i6,a,i6,a,i8)') 344 . 'superc: Internal auxiliary supercell:', 345 . nsc(1), ' x', nsc(2), ' x', nsc(3), ' =', ncells 346 347 write(6,'(a,1x,i5,2(1x,i6))') 348 $ 'superc: Number of atoms, orbitals, and projectors: ', 349 $ na_s, no_s, nokb_s 350 351 endif 352 353 end subroutine superc 354 355 SUBROUTINE SUPERX( UCELL, NSC, NA, MAXA, XA, SCELL ) 356 357C ********************************************************************** 358C Generates supercell vectors and atomic positions. 359C Written by J.M.Soler, August 1998 360C *************** Input ************************************************ 361C Real*8 UCELL(3,3) : Unit cell vectors UCELL(Ixyz,Ivector) 362C Integer NSC(3) : Number of cells in each supercell direction: 363C SCELL(ix,i) = UCELL(ix,i) * NSC(i) 364C Integer NA : Number of atoms in unit cell 365C Integer MAXA : Second dimension of XA 366C *************** Input and output ************************************* 367C Real*8 XA(3,MAXA) : Atomic positions in unit cell (input) and 368C in supercell (output), in cartesian coord. 369C Real*8 SCELL(3,3) : Supercell vectors 370C *********** Units **************************************************** 371C Units of CELL and XA are arbitrary but must be the same 372C *********** Behavior ************************************************* 373C - If NA*NCELLS > MAXA (where NCELLS is the total number of cells), 374C the supercell atomic coordinates are not generated. 375C - The first supercell atoms are those of the initial unit cell, i.e. 376C the positions XA(i,ia) for (ia.le.NA) are not modified. 377C - The remaining atoms are ordered by unit cells, i.e. the atom ia 378C is equivalent to the unit-cell atom ja=MOD(ia-1,NA)+1 379C ********************************************************************** 380 381 IMPLICIT NONE 382 INTEGER MAXA, NA, NSC(3) 383 DOUBLE PRECISION SCELL(3,3), UCELL(3,3), XA(3,MAXA) 384 385C Internal variables 386 INTEGER I, I1, I2, I3, IA, IX, JA, NCELLS 387 DOUBLE PRECISION XC(3) 388 389C Find supercell vectors 390 DO 10 I = 1,3 391 DO 5 IX = 1,3 392 SCELL(IX,I) = UCELL(IX,I) * NSC(I) 393 5 CONTINUE 394 10 CONTINUE 395 396C Expand atomic positions to supercell 397 NCELLS = NSC(1) * NSC(2) * NSC(3) 398 IF (NA*NCELLS .LE. MAXA) THEN 399 IA = 0 400 DO 60 I3 = 0,NSC(3)-1 401 DO 50 I2 = 0,NSC(2)-1 402 DO 40 I1 = 0,NSC(1)-1 403 DO 15 IX = 1,3 404 XC(IX) = UCELL(IX,1)*I1 + UCELL(IX,2)*I2 + UCELL(IX,3)*I3 405 15 CONTINUE 406 DO 30 JA = 1,NA 407 IA = IA + 1 408 DO 20 IX = 1,3 409 XA(IX,IA) = XA(IX,JA) + XC(IX) 410 20 CONTINUE 411 30 CONTINUE 412 40 CONTINUE 413 50 CONTINUE 414 60 CONTINUE 415 ENDIF 416 417 END subroutine superx 418 419 end module atomlist 420 421 422 423