1 Subroutine grid_atom_type_info 2c 3c$Id$ 4c 5 implicit none 6#include "errquit.fh" 7c 8#include "inp.fh" 9#include "util.fh" 10#include "global.fh" 11#include "stdio.fh" 12#include "cdft.fh" 13#include "geom.fh" 14#include "mafdecls.fh" 15c 16 double precision BSrad(105) 17 double precision ictr_coord(3), ictr_chg 18 double precision autoang, angtoau, EPS 19 parameter (autoang = 0.529177249d0, angtoau = 1.d0/autoang) 20 parameter (EPS = 1.d-20) 21 integer itype, ictr, iaz, i_atomic_number, icenter, jcenter 22c 23 integer lcoord, icoord, lcharge, icharge, ltags, itags,iptr 24 logical same_atom, same_bq, isbq 25c 26 character*16 element 27 character*16 tag 28 character*2 symbol 29c 30 logical lnewtype 31 logical atom_tag_check 32 logical oprint, oprint_grid 33c 34 external atom_tag_check 35c 36c 37c Table of Bragg-Slater Atomic Radii (Angstroms) 38c 39c Bragg-Slater radii: J.C. Slater, Symmetry and Energy Bands in Crystals, 40c Dover, N.Y. 1972, page 55. 41c The radii of noble gas atoms are set to be equal 42c to the radii of the corresponding halogen atoms. 43c The radius of At is set to be equal to the radius of 44c Po. 45c The radius of Fr is set to be equal to the radius of 46c Cs. 47c 48c H He Li Be B C N O F Ne 49 Data BSrad/0.35,0.35,1.45,1.05,0.85,0.70,0.65,0.60,0.50,0.50, 50c Na Mg Al Si P S Cl Ar K Ca 51 & 1.80,1.50,1.25,1.10,1.00,1.00,1.00,1.00,2.20,1.80, 52c Sc Ti V Cr Mn Fe Co Ni Cu Zn 53 & 1.60,1.40,1.35,1.40,1.40,1.40,1.35,1.35,1.35,1.35, 54c Ga Ge As Se Br Kr Rb Sr Y Zr 55 & 1.30,1.25,1.15,1.15,1.15,1.15,2.35,2.00,1.80,1.55, 56c Nb Mo Tc Ru Rh Pd Ag Cd In Sn 57 & 1.45,1.45,1.35,1.30,1.35,1.40,1.60,1.55,1.55,1.45, 58c Sb Te I Xe Cs Ba La Ce Pr Nd 59 & 1.45,1.40,1.40,1.40,2.60,2.15,1.95,1.85,1.85,1.85, 60c Pm Sm Eu Gd Tb Dy Ho Er Tm Yb 61 & 1.85,1.85,1.85,1.80,1.75,1.75,1.75,1.75,1.75,1.75, 62c Lu Hf Ta W Re Os Ir Pt Au Hg 63 & 1.75,1.55,1.45,1.35,1.35,1.30,1.35,1.35,1.35,1.50, 64c Tl Pb Bi Po At Rn Fr Ra Ac Th 65 & 1.90,1.80,1.60,1.90,1.90,1.90,2.60,2.15,1.95,1.80, 66c Pa U Np Pu Am Cm Bk Cf Es Fm 67 & 1.80,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75, 68c Md No Lr Unq Unp 69 & 1.75,1.75,1.75,1.55,1.55/ 70c 71c Set print options. 72c 73 oprint = util_print('quadrature', print_high) 74 oprint_grid = util_print('griddebug', print_debug) 75c 76c allocate space for atomic coordinates and charges 77c 78 if (.not. Ma_Push_Get(MT_Dbl,ncenters*3,'coordinates',lcoord, 79 & icoord))call errquit( 80 G 'grid_atom_type_info: failed to alloc coordinates',0, MA_ERR) 81 if (.not. Ma_Push_Get(MT_Dbl,ncenters,'charges',lcharge, 82 & icharge))call errquit( 83 G 'grid_atom_type_info: failed to alloc charges',0, MA_ERR) 84 if (.not. Ma_Push_Get(MT_Byte, ncenters*16, 'center tags', 85 & ltags, itags))call errquit( 86 G 'grid_atom_type_info: failed to alloc center tags',0, 87 M MA_ERR) 88c 89c Get ncenter tags, coordinates, and charges from the geometry object. 90c 91 if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags), 92 & Dbl_MB(icoord), Dbl_MB(icharge))) 93 & call errquit('gridatom_type_info: geom_cart_get failed',geom, 94 & GEOM_ERR) 95c 96c generate number of atom types and atom type array iatype(icenter) 97c 98 ntypes = 0 99 do icenter = 1, ncenters 100c 101c is this a new type of atom? 102c 103 lnewtype = .true. 104 isbq = geom_isbq(geom,icenter) 105 do jcenter = 1, icenter - 1 106 same_atom = Dbl_MB(icharge + icenter - 1) .eq. 107 & Dbl_MB(icharge + jcenter - 1) 108 same_bq = geom_isbq(geom,jcenter) .and. isbq 109 same_atom = same_atom .or. same_bq 110 if (same_atom .and. 111 & atom_tag_check(Byte_MB(itags + (icenter - 1)*16), 112 & Byte_MB(itags + (jcenter - 1)*16)) 113 & )then ! same atom type 114 lnewtype = .false. 115 iatype(icenter) = iatype(jcenter) 116 goto 100 117 endif 118 enddo 119 100 continue 120 if (lnewtype)then 121 ntypes = ntypes + 1 122 iatype(icenter) = ntypes 123 endif 124 enddo 125 if (ntypes.gt.dft_ntags_bsmx)then 126 write(LuOut,*) 'grid_atom_type_info: Too many types of atoms.' 127 call errquit(' grid_atom_type_info: raise dft_ntags_bsmx',2, 128 & GEOM_ERR) 129 end if 130c 131c set up type-indexed znuc array; znuc_atom_type 132c 133 do itype = 1, ntypes 134 do icenter = 1, ncenters 135 if (iatype(icenter) .eq. itype)then 136c 137c center icenter is of type itype; assign charge 138c 139 znuc_atom_type(itype) = Dbl_MB(icharge + icenter - 1) 140 goto 110 ! next type 141 endif 142 enddo 143 110 continue 144 enddo 145c 146c Define the atomic Bragg-Slater radii for each atom type. 147c 148 do 50 itype = 1, ntypes 149c 150c find an atom of this kind in the complete list 151c 152 do ictr = 1, ncenters 153 if (iatype(ictr).eq.itype) then 154 iaz = ictr 155 if (.not. geom_cent_get(geom, ictr, tag, 156 & ictr_coord, ictr_chg))call errquit 157 & ('grid_atom_type_info: geom_cent_get failed', 0, 158 & GEOM_ERR) 159 goto 40 160 endif 161 enddo 162 40 continue 163c 164 if (abs(znuc_atom_type(itype)).lt.EPS) then ! uncharged ghost atom; 165c 166c identify atom label following "bq" or "x" (we already know this 167c cannot be element Xeon). 168c 169 iptr=3 170c hack for nbo 171 if(tag(3:4).eq.'gh') iptr=5 172c 173c Dummy centers (i.e. X, XH, X1, whatever) by convention 174c cannot have charges or basis sets. They are purely 175c mathematical constructs to specify complex Z-matrices. 176c Hence they should always have i_atomic_number .eq. 0 177c no matter what follows the X. 178c 179 if(tag(1:1).eq.'X'.or.tag(1:1).eq.'x') iptr=2 180 if (.not. geom_tag_to_element(tag(iptr:), symbol, 181 & element, i_atomic_number)) then 182 if (inp_compare(.false.,tag(1:2),'bq')) then 183 if(bqdontcare) then 184 i_atomic_number = 1 185 else 186 i_atomic_number = 0 187 endif 188 elseif (inp_compare(.false.,tag(1:1),'X')) then 189 i_atomic_number = 0 190 else 191 call errquit( 192 & 'grid_atom_type_info: non-bq center with zero charge', 193 & 0, INPUT_ERR) 194 & 195 endif 196 else 197 if (inp_compare(.false.,tag(1:1),'X')) then 198 i_atomic_number = 0 199 endif 200 endif 201c 202c 203 if (i_atomic_number.eq.0)then 204 bsrad_atom_type(itype) = EPS 205 else 206 bsrad_atom_type(itype) = BSrad(i_atomic_number)*ANGTOAU 207 endif 208 else ! center is charged 209c 210c no quadrature grids on charged ghost atoms 211c 212 if (.not. geom_tag_to_element(tag, symbol, 213 & element, i_atomic_number)) then 214 if (symbol .ne. 'bq') call errquit 215 & ('grid_atom_type_info: center is neither atom nor bq', 216 & 0, INPUT_ERR) 217 endif 218c 219 if (i_atomic_number.ne.0)then ! not ghost atom 220 ityp2ctr(itype)=ictr 221c 222 bsrad_atom_type(itype) = BSrad(i_atomic_number)*ANGTOAU 223 if (bsrad_atom_type(itype).lt.EPS)then ! no radius found for atom 224 write(LuOut,*)' index ', 225 & int(abs(znuc_atom_type(itype)) + EPS) 226 write(LuOut,*)' BSR ', 227 & BSrad(int(abs(znuc_atom_type(itype)) + EPS)) 228 write(LuOut,*)' grid_atom_type_info: ', 229 & ' Undefined atomic radius ' 230 write(LuOut,*)' for atom type', itype 231 call errquit('Exiting in grid_atom_type_info.',1, 232 & UNKNOWN_ERR) 233 endif 234 else ! atomic number zero; charged ghost atom 235 bsrad_atom_type(itype) = EPS 236 endif 237 endif 238 50 continue 239c 240c Build logical vector to tag centers as point charge centers or 241c real centers with basis functions 242c 243 do 55 ictr = 1, ncenters 244 itype = iatype(ictr) 245 if (bsrad_atom_type(itype).le.EPS)then 246 iatype_pt_chg(ictr) = .true. 247 else 248 iatype_pt_chg(ictr) = .false. 249 endif 250 55 continue 251 if (.not. MA_Pop_Stack(ltags)) 252 & call errquit('grid_atom_type_info: pop stack failed.',0, 253 & MA_ERR) 254 if (.not. MA_Pop_Stack(lcharge)) 255 & call errquit('grid_atom_type_info: pop stack failed.',0, 256 & MA_ERR) 257 if (.not. MA_Pop_Stack(lcoord)) 258 & call errquit('grid_atom_type_info: pop stack failed.',0, 259 & MA_ERR) 260c 261c debug writes 262c 263 if (ga_nodeid().eq.0.and.oprint_grid)then 264 write(LuOut,*)' iatype(ncenters) ', 265 & (iatype(ictr),ictr = 1, ncenters) 266 write(LuOut,*)' iatype_pt_chg(ncenters) ', 267 & (iatype_pt_chg(ictr),ictr = 1, ncenters) 268 write(LuOut,*)' bsrad_atom_type(ntypes) ', 269 & (bsrad_atom_type(itype),itype = 1, ntypes) 270 write(LuOut,*)' znuc_atom_type(ntypes) ', 271 & (znuc_atom_type(itype),itype = 1, ntypes) 272 endif 273 return 274 end 275 276 logical function atom_tag_check(atom_tag_i, atom_tag_j) 277c 278 implicit none 279c 280 character*16 atom_tag_i, atom_tag_j 281c write(*,*)' atom_tag_i = ', atom_tag_i 282c write(*,*)' atom_tag_j = ', atom_tag_j 283 if (atom_tag_i .eq. atom_tag_j)then 284 atom_tag_check = .true. 285 else 286 atom_tag_check = .false. 287 endif 288 return 289 end 290 291cc AJL/Begin/FDE 292c 293 Subroutine grid_atom_type_info_fde 294c 295c$Id$ 296c 297 implicit none 298#include "errquit.fh" 299c 300#include "inp.fh" 301#include "util.fh" 302#include "global.fh" 303#include "stdio.fh" 304#include "cdft.fh" 305#include "geom.fh" 306#include "mafdecls.fh" 307c 308 double precision BSrad(105) 309 double precision ictr_coord(3), ictr_chg 310 double precision autoang, angtoau, EPS 311 parameter (autoang = 0.529177249d0, angtoau = 1.d0/autoang) 312 parameter (EPS = 1.d-20) 313 integer itype, ictr, iaz, i_atomic_number, icenter, jcenter 314c 315 integer lcoord, icoord, lcharge, icharge, ltags, itags,iptr 316 logical same_atom, same_bq, isbq 317c 318 character*16 element 319 character*16 tag 320 character*2 symbol 321c 322 logical lnewtype 323 logical atom_tag_check 324 logical oprint, oprint_grid 325c AJL/Begin 326 integer ncenters_fde 327 integer ntypes_fde 328c AJL/End 329c 330 External atom_tag_check 331c 332c 333c Table of Bragg-Slater Atomic Radii (Angstroms) 334c 335c Bragg-Slater radii: J.C. Slater, Symmetry and Energy Bands in 336Crystals, 337c Dover, N.Y. 1972, page 55. 338c The radii of noble gas atoms are set to be equal 339c to the radii of the corresponding halogen atoms. 340c The radius of At is set to be equal to the radius of 341c Po. 342c The radius of Fr is set to be equal to the radius of 343c Cs. 344c 345c H He Li Be B C N O F Ne 346 Data BSrad/0.35,0.35,1.45,1.05,0.85,0.70,0.65,0.60,0.50,0.50, 347c Na Mg Al Si P S Cl Ar K Ca 348 & 1.80,1.50,1.25,1.10,1.00,1.00,1.00,1.00,2.20,1.80, 349c Sc Ti V Cr Mn Fe Co Ni Cu Zn 350 & 1.60,1.40,1.35,1.40,1.40,1.40,1.35,1.35,1.35,1.35, 351c Ga Ge As Se Br Kr Rb Sr Y Zr 352 & 1.30,1.25,1.15,1.15,1.15,1.15,2.35,2.00,1.80,1.55, 353c Nb Mo Tc Ru Rh Pd Ag Cd In Sn 354 & 1.45,1.45,1.35,1.30,1.35,1.40,1.60,1.55,1.55,1.45, 355c Sb Te I Xe Cs Ba La Ce Pr Nd 356 & 1.45,1.40,1.40,1.40,2.60,2.15,1.95,1.85,1.85,1.85, 357c Pm Sm Eu Gd Tb Dy Ho Er Tm Yb 358 & 1.85,1.85,1.85,1.80,1.75,1.75,1.75,1.75,1.75,1.75, 359c Lu Hf Ta W Re Os Ir Pt Au Hg 360 & 1.75,1.55,1.45,1.35,1.35,1.30,1.35,1.35,1.35,1.50, 361c Tl Pb Bi Po At Rn Fr Ra Ac Th 362 & 1.90,1.80,1.60,1.90,1.90,1.90,2.60,2.15,1.95,1.80, 363c Pa U Np Pu Am Cm Bk Cf Es Fm 364 & 1.80,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75,1.75, 365c Md No Lr Unq Unp 366 & 1.75,1.75,1.75,1.55,1.55/ 367c 368c Set print options. 369c 370 oprint = util_print('quadrature', print_high) 371 oprint_grid = util_print('griddebug', print_debug) 372c 373c Get the number of FDE centers 374c 375 if (.not. geom_ncent(geom_fde, ncenters_fde)) 376 & call errquit('grid_atom_type_info_fde: 377 & geom_ncent failed',73, GEOM_ERR) 378c 379c allocate space for atomic coordinates and charges 380c 381 if (.not. Ma_Push_Get(MT_Dbl,ncenters_fde*3,'coordinates',lcoord, 382 & icoord))call errquit( 383 . 'grid_atom_type_info: failed to alloc coordinates',0, MA_ERR) 384 if (.not. Ma_Push_Get(MT_Dbl,ncenters_fde,'charges',lcharge, 385 & icharge))call errquit( 386 ' 'grid_atom_type_info: failed to alloc charges',0, MA_ERR) 387 if (.not. Ma_Push_Get(MT_Byte, ncenters_fde*16, 'center tags', 388 & ltags, itags))call errquit( 389 . 'grid_atom_type_info: failed to alloc center tags',0, MA_ERR) 390c 391c Get ncenter tags, coordinates, and charges from the geometry object. 392c 393 if (.not. geom_cart_get(geom_fde, ncenters_fde, Byte_MB(itags), 394 & Dbl_MB(icoord), Dbl_MB(icharge))) 395 & call errquit('gridatom_type_info: geom_cart_get failed',geom_fde, 396 & GEOM_ERR) 397c 398c generate number of atom types and atom type array iatype(icenter) 399c 400 ntypes_fde = 0 401 do icenter = 1, ncenters_fde 402c 403c is this a new type of atom? 404c 405 lnewtype = .true. 406 isbq = geom_isbq(geom_fde,icenter) 407 do jcenter = 1, icenter - 1 408 same_atom = Dbl_MB(icharge + icenter - 1) .eq. 409 & Dbl_MB(icharge + jcenter - 1) 410 same_bq = geom_isbq(geom_fde,jcenter) .and. isbq 411 same_atom = same_atom .or. same_bq 412 if (same_atom .and. 413 & atom_tag_check(Byte_MB(itags + (icenter - 1)*16), 414 & Byte_MB(itags + (jcenter - 1)*16)) 415 & )then ! same atom type 416 lnewtype = .false. 417 iatype_fde(icenter) = iatype_fde(jcenter) 418 goto 100 419 endif 420 enddo 421 100 continue 422 if (lnewtype)then 423 ntypes_fde = ntypes_fde + 1 424 iatype_fde(icenter) = ntypes_fde 425 endif 426 enddo 427 if (ntypes_fde.gt.dft_ntags_bsmx)then 428 write(LuOut,*) 'grid_atom_type_info_fde: Too many types of atoms.' 429 call errquit(' grid_atom_type_info_fde: raise dft_ntags_bsmx',2, 430 & GEOM_ERR) 431 end if 432c 433c set up type-indexed znuc array; znuc_atom_type 434c 435 do itype = 1, ntypes_fde 436 do icenter = 1, ncenters_fde 437 if (iatype_fde(icenter) .eq. itype)then 438c 439c center icenter is of type itype; assign charge 440c 441 znuc_atom_type_fde(itype) = Dbl_MB(icharge + icenter - 1) 442 goto 110 ! next type 443 endif 444 enddo 445 110 continue 446 enddo 447c 448c Define the atomic Bragg-Slater radii for each atom type. 449c 450 do 50 itype = 1, ntypes_fde 451c 452c find an atom of this kind in the complete list 453c 454 do ictr = 1, ncenters_fde 455 if (iatype_fde(ictr).eq.itype) then 456 iaz = ictr 457 if (.not. geom_cent_get(geom_fde, ictr, tag, 458 & ictr_coord, ictr_chg))call errquit 459 & ('grid_atom_type_info_fde: geom_cent_get failed', 0, 460 & GEOM_ERR) 461 goto 40 462 endif 463 enddo 464 40 continue 465c 466 if (abs(znuc_atom_type_fde(itype)).lt.EPS) then ! uncharged ghost atom; 467c 468c identify atom label following "bq" or "x" (we already know this 469c cannot be element Xeon). 470c 471 iptr=3 472c hack for nbo 473 if(tag(3:4).eq.'gh') iptr=5 474c 475c Dummy centers (i.e. X, XH, X1, whatever) by convention 476c cannot have charges or basis sets. They are purely 477c mathematical constructs to specify complex Z-matrices. 478c Hence they should always have i_atomic_number .eq. 0 479c no matter what follows the X. 480c 481 if(tag(1:1).eq.'X'.or.tag(1:1).eq.'x') iptr=2 482 if (.not. geom_tag_to_element(tag(iptr:), symbol, 483 & element, i_atomic_number)) then 484 if (inp_compare(.false.,tag(1:2),'bq')) then 485 if(bqdontcare) then 486 i_atomic_number = 1 487 else 488 i_atomic_number = 0 489 endif 490 elseif (inp_compare(.false.,tag(1:1),'X')) then 491 i_atomic_number = 0 492 else 493 call errquit('grid_atom_type_info_fde: 494 & non-bq center with zero charge', 0, INPUT_ERR) 495 endif 496 else 497 if (inp_compare(.false.,tag(1:1),'X')) then 498 i_atomic_number = 0 499 endif 500 endif 501c 502c 503 if (i_atomic_number.eq.0)then 504 bsrad_atom_type_fde(itype) = EPS 505 else 506 bsrad_atom_type_fde(itype) = 507 & BSrad(i_atomic_number)*ANGTOAU 508 endif 509 else ! center is charged 510c 511c no quadrature grids on charged ghost atoms 512c 513 if (.not. geom_tag_to_element(tag, symbol, 514 & element, i_atomic_number)) then 515 if (symbol .ne. 'bq') call errquit 516 & ('grid_atom_type_info_fde: center is neither atom nor 517 & bq', 0, INPUT_ERR) 518 endif 519c 520 if (i_atomic_number.ne.0)then ! not ghost atom 521c ityp2ctr(itype)=ictr 522c 523 bsrad_atom_type_fde(itype) = 524 & BSrad(i_atomic_number)*ANGTOAU 525 if (bsrad_atom_type_fde(itype).lt.EPS)then ! no radius found for atom 526 write(LuOut,*)' index ', 527 & int(abs(znuc_atom_type_fde(itype)) + EPS) 528 write(LuOut,*)' BSR ', 529 & BSrad(int(abs(znuc_atom_type_fde(itype)) + EPS)) 530 write(LuOut,*)' grid_atom_type_info_fde: ', 531 & ' Undefined atomic radius ' 532 write(LuOut,*)' for atom type', itype 533 call errquit('Exiting in grid_atom_type_info_fde.',1, 534 & UNKNOWN_ERR) 535 endif 536 else ! atomic number zero; charged ghost atom 537 bsrad_atom_type_fde(itype) = EPS 538 endif 539 endif 540 50 continue 541c 542c Build logical vector to tag centers as point charge centers or 543c real centers with basis functions 544c 545 do 55 ictr = 1, ncenters_fde 546 itype = iatype_fde(ictr) 547 if (bsrad_atom_type_fde(itype).le.EPS)then 548 iatype_pt_chg_fde(ictr) = .true. 549 else 550 iatype_pt_chg_fde(ictr) = .false. 551 endif 552 55 continue 553 554 if (.not. MA_Pop_Stack(ltags)) 555 & call errquit('grid_atom_type_info_fde: pop stack failed.',0, 556 & MA_ERR) 557 if (.not. MA_Pop_Stack(lcharge)) 558 & call errquit('grid_atom_type_info_fde: pop stack failed.',0, 559 & MA_ERR) 560 if (.not. MA_Pop_Stack(lcoord)) 561 & call errquit('grid_atom_type_info_fde: pop stack failed.',0, 562 & MA_ERR) 563c 564c debug writes 565c 566 if (ga_nodeid().eq.0.and.oprint_grid)then 567 write(LuOut,*)' iatype_fde(ncenters_fde) ', 568 & (iatype(ictr),ictr = 1, ncenters_fde) 569 write(LuOut,*)' iatype_pt_chg_fde(ncenters_fde) ', 570 & (iatype_pt_chg(ictr),ictr = 1, ncenters_fde) 571 write(LuOut,*)' bsrad_atom_type_fde(ntypes_fde) ', 572 & (bsrad_atom_type(itype),itype = 1, ntypes_fde) 573 write(LuOut,*)' znuc_atom_type_fde(ntypes_fde) ', 574 & (znuc_atom_type(itype),itype = 1, ntypes_fde) 575 endif 576 return 577 end 578c 579cc AJL/END 580