1 2* 3* $Id$ 4* 5 6* *********************************** 7* * * 8* * borbs_init * 9* * * 10* *********************************** 11 logical function borbs_init() 12 implicit none 13 14#include "bafdecls.fh" 15#include "borbs.fh" 16 17* **** local variables **** 18 integer npack1,nion 19 logical value 20 21* **** external functions ***** 22 integer ion_nkatm_qm,ion_nion_qm 23 external ion_nkatm_qm,ion_nion_qm 24 25 26 call Pack_npack(1,npack1) 27 nion = ion_nion_qm() 28 nkatmx = ion_nkatm_qm() 29 30* **** allocate borb datastructure and borb tag lists**** 31 call borb_projector_init(2*nkatmx) 32 value = BA_alloc_get(mt_int,(nkatmx),'borbs',borbs(2),borbs(1)) 33 34 value = value.and. 35 > BA_alloc_get(mt_int,(nkatmx),'lmmax',lmmax(2),lmmax(1)) 36 value = value.and. 37 > BA_alloc_get(mt_int,(nkatmx),'lmax',lmax(2),lmax(1)) 38 value = value.and. 39 > BA_alloc_get(mt_int,(nkatmx),'locp',locp(2),locp(1)) 40 value = value.and. 41 > BA_alloc_get(mt_dbl,(nkatmx),'rcut',rcut(2),rcut(1)) 42 value = value.and. 43 > BA_alloc_get(mt_dbl,(nkatmx),'lmbda',lmbda(2),lmbda(1)) 44 45 value = value.and. 46 > BA_alloc_get(mt_int,(nion*norbs_max), 47 > 'lmborb',lmborb(2),lmborb(1)) 48 value = value.and. 49 > BA_alloc_get(mt_int,(nion*norbs_max), 50 > 'iaborb',iaborb(2),iaborb(1)) 51 value = value.and. 52 > BA_alloc_get(mt_int,(nion*norbs_max), 53 > 'iiborb',iiborb(2),iiborb(1)) 54 value = value.and. 55 > BA_alloc_get(mt_int,(nion*norbs_max), 56 > 'basisborb',basisborb(2),basisborb(1)) 57 58 borbs_init = value 59 return 60 end 61 62 63* *********************************** 64* * * 65* * borbs_end * 66* * * 67* *********************************** 68 subroutine borbs_end() 69 implicit none 70 71#include "bafdecls.fh" 72#include "errquit.fh" 73#include "borbs.fh" 74 75* **** local variables **** 76 logical value 77 78* **** deallocate projector data **** 79 call borb_projector_end() 80 81 value = BA_free_heap(borbs(2)) 82 value = value.and.BA_free_heap(lmmax(2)) 83 value = value.and.BA_free_heap(lmax(2)) 84 value = value.and.BA_free_heap(locp(2)) 85 value = value.and.BA_free_heap(rcut(2)) 86 value = value.and.BA_free_heap(lmbda(2)) 87 value = value.and.BA_free_heap(lmborb(2)) 88 value = value.and.BA_free_heap(iaborb(2)) 89 value = value.and.BA_free_heap(iiborb(2)) 90 value = value.and.BA_free_heap(basisborb(2)) 91 if (.not. value) call errquit('borbs_end:freeing heap memory',0, 92 & MA_ERR) 93 94 return 95 end 96 97 98* *********************************** 99* * * 100* * borbs_norbs * 101* * * 102* *********************************** 103 104 integer function borbs_norbs(ia) 105 implicit none 106 integer ia 107 108#include "bafdecls.fh" 109#include "borbs.fh" 110 111 borbs_norbs = int_mb(lmmax(1)+ia-1) 112 return 113 end 114 115 116* *********************************** 117* * * 118* * borbs_nbasis * 119* * * 120* *********************************** 121 integer function borbs_nbasis() 122 implicit none 123 124#include "borbs.fh" 125 126 borbs_nbasis = ibasis 127 return 128 end 129 130* *********************************** 131* * * 132* * borbs_lmax * 133* * * 134* *********************************** 135 136 integer function borbs_lmax(ia) 137 implicit none 138 integer ia 139 140#include "bafdecls.fh" 141#include "borbs.fh" 142 143 borbs_lmax = int_mb(lmax(1)+ia-1)-1 144 return 145 end 146 147 148* *********************************** 149* * * 150* * borbs_rcut * 151* * * 152* *********************************** 153 real*8 function borbs_rcut(ia) 154 implicit none 155 integer ia 156 157#include "bafdecls.fh" 158#include "borbs.fh" 159 160 borbs_rcut = dbl_mb(rcut(1)+ia-1) 161 return 162 end 163 164* *********************************** 165* * * 166* * borbs_lmbda * 167* * * 168* *********************************** 169 real*8 function borbs_lmbda(ia) 170 implicit none 171 integer ia 172 173#include "bafdecls.fh" 174#include "borbs.fh" 175 176 borbs_lmbda = dbl_mb(lmbda(1)+ia-1) 177 return 178 end 179 180 181* *********************************** 182* * * 183* * borbs_l * 184* * * 185* *********************************** 186 187 integer function borbs_l(ia,n) 188 implicit none 189 integer ia 190 integer n ! basis number 191 192#include "bafdecls.fh" 193#include "borbs.fh" 194 195* *** local variables *** 196 integer l,m,lm 197 198 lm = int_mb(lmborb(1)+n-1) 199 l = 0 200 if (lm.eq.1) l = 0 201 202 if (lm.eq.2) l = 1 203 if (lm.eq.3) l = 1 204 if (lm.eq.4) l = 1 205 206 if (lm.eq.5) l = 2 207 if (lm.eq.6) l = 2 208 if (lm.eq.7) l = 2 209 if (lm.eq.8) l = 2 210 if (lm.eq.9) l = 2 211 212 if (lm.eq.10) l = 3 213 if (lm.eq.11) l = 3 214 if (lm.eq.12) l = 3 215 if (lm.eq.13) l = 3 216 if (lm.eq.14) l = 3 217 if (lm.eq.15) l = 3 218 if (lm.eq.16) l = 3 219 220 borbs_l = l 221 return 222 end 223 224 225* *********************************** 226* * * 227* * borbs_get_basis_number * 228* * * 229* *********************************** 230 231 integer function borbs_get_basis_number(ii,lm) 232 implicit none 233 integer ii,lm 234 235#include "bafdecls.fh" 236#include "borbs.fh" 237 238 borbs_get_basis_number=int_mb(basisborb(1)+lm-1+(ii-1)*norbs_max) 239 240 return 241 end 242 243 244 245* *********************************** 246* * * 247* * borbs_normalize * 248* * * 249* *********************************** 250 subroutine borbs_normalize() 251 implicit none 252 253#include "bafdecls.fh" 254#include "borbs.fh" 255 256* **** local variables **** 257 integer lm,ia,npack1,nbrillq,nbq,shift 258 real*8 sum 259 260* **** external functions **** 261 integer Pneb_nbrillq,borb_projector_get_ptr 262 external Pneb_nbrillq,borb_projector_get_ptr 263 264 265 nbrillq = Pneb_nbrillq() 266 267* **** Normalize atomic orbitals in k space ***** 268 do nbq=1,nbrillq 269 do ia=1,nkatmx 270 do lm=1,int_mb(lmmax(1)+ia-1) 271 272 shift = borb_projector_get_ptr(int_mb(borbs(1)+ia-1),nbq,lm) 273 call Cram_rr_dot(nbq,dbl_mb(shift),dbl_mb(shift),sum) 274 sum = 1.0d0/dsqrt(sum) 275 call Cram_r_SMul1(nbq,sum,dbl_mb(shift)) 276 end do 277 end do 278 end do 279 280 return 281 end 282 283* *********************************** 284* * * 285* * borbs_weight * 286* * * 287* *********************************** 288 real*8 function borbs_weight(n) 289 implicit none 290 integer n ! basis number 291 292#include "bafdecls.fh" 293#include "borbs.fh" 294 295* **** local variables **** 296 integer ia 297 real*8 zv,zcount 298 299* **** external functions **** 300 real*8 psp_zv 301 external psp_zv 302 303 ia = int_mb(iaborb(1)+n-1) 304 zcount = int_mb(lmmax(1)+ia-1) 305 zv = psp_zv(ia) 306 307 borbs_weight = (zv/zcount) 308 return 309 end 310 311* *********************************** 312* * * 313* * borbs_borb * 314* * * 315* *********************************** 316 317 subroutine borbs_borb(nbq,n,borb) 318 implicit none 319 integer nbq 320 integer n ! basis number 321 complex*16 borb(*) ! return orbital 322 323#include "bafdecls.fh" 324#include "errquit.fh" 325#include "borbs.fh" 326 327* **** local variables **** 328 logical value 329 integer lm,ia,ii 330 integer nfft3d,npack1,shift 331 integer exi(2) 332 complex*16 cxr 333 334* **** external functions **** 335 logical is_sORd 336 external is_sORd 337 integer borb_projector_get_ptr 338 external borb_projector_get_ptr 339 340 call C3dB_nfft3d(1,nfft3d) 341 call Cram_npack(nbq,npack1) 342 343 value = BA_push_get(mt_dcpl,npack1,'exi', exi(2), exi(1)) 344 if (.not. value) call errquit('borbs_borb:out of heap memory',0, 345 & MA_ERR) 346 347* **** structure factor **** 348 lm = int_mb(lmborb(1)+n-1) 349 ia = int_mb(iaborb(1)+n-1) 350 ii = int_mb(iiborb(1)+n-1) 351 call cstrfac_pack(nbq,ii,dcpl_mb(exi(1))) 352 call cstrfac_k(ii,nbq,cxr) 353 call zscal(npack1,cxr,dcpl_mb(exi(1)),1) 354 355 shift = borb_projector_get_ptr(int_mb(borbs(1)+ia-1),nbq,lm) 356 357 358* **** phase factor does not matter therefore **** 359* **** (-i)^l is the same as (i)^l in the **** 360* **** Rayleigh scattering formula **** 361 362* *** current function is s or d **** 363 if (is_sORd(lm,int_mb(lmax(1)+ia-1), 364 > int_mb(locp(1)+ia-1)) 365 > ) then 366 call Cram_rc_Mul(nbq,dbl_mb(shift),dcpl_mb(exi(1)),borb) 367 368* *** current function is p or f **** 369 else 370 call Cram_irc_Mul(nbq,dbl_mb(shift),dcpl_mb(exi(1)),borb) 371 end if 372 373 374 value = BA_pop_stack(exi(2)) 375 if (.not. value) call errquit('borbs_borb:freeing heap memory',0, 376 & MA_ERR) 377 378 return 379 end 380 381 382* *********************************** 383* * * 384* * borbs_read * 385* * * 386* *********************************** 387 subroutine borbs_read(fname, 388 > version, 389 > nfft,unita, 390 > atom, 391 > lmmax,lmax,locp,rcut,lmbda, 392 > npack1,borbs_tag, 393 > tmp,tmp2, 394 > ierr) 395 implicit none 396 character*50 fname 397 integer version 398 integer nfft(3) 399 real*8 unita(3,3) 400 character*2 atom 401 integer lmmax,lmax,locp 402 real*8 rcut,lmbda 403 integer npack1 404 integer borbs_tag 405 complex*16 tmp(*) 406 real*8 tmp2(*) 407 integer ierr 408 409#include "bafdecls.fh" 410#include "btdb.fh" 411#include "util.fh" 412 413* *** local variables *** 414 logical mprint,value 415 integer MASTER,taskid,taskid_k 416 parameter(MASTER=0) 417 integer n,l,nbrillioun,nb,nbq,pk 418 integer msglen 419 integer iatom(2) 420 character*255 full_filename 421 real*8 kv(3) 422 423 424* ***** external functions **** 425 integer brillioun_nbrillioun,borb_projector_alloc 426 external brillioun_nbrillioun,borb_projector_alloc 427 integer brillioun_nbrillq 428 external brillioun_nbrillq 429 real*8 brillioun_all_k 430 external brillioun_all_k 431 logical control_print 432 external control_print 433 434 call Parallel_taskid(taskid) 435 call Parallel3d_taskid_k(taskid_k) 436 mprint = (taskid.eq.MASTER).and.control_print(print_medium) 437 438 439* **** open fname binary file **** 440 if (taskid.eq.MASTER) then 441 call util_file_name_noprefix(fname,.false., 442 > .false., 443 > full_filename) 444 l = index(full_filename,' ') - 1 445 call openfile(5,full_filename,l,'r',l) 446 call iread(5,version,1) 447 call iread(5,nfft,3) 448 call dread(5,unita,9) 449 call cread(5,atom,2) 450 call iread(5,lmax,1) 451 call iread(5,locp,1) 452 call dread(5,rcut,1) 453 call dread(5,lmbda,1) 454 455 call iread(5,nbrillioun,1) 456 ierr = 0 457 if (nbrillioun.eq.brillioun_nbrillioun()) then 458 do nb=1,nbrillioun 459 call dread(5,kv,3) 460 if ((brillioun_all_k(1,nb).ne.kv(1)).or. 461 > (brillioun_all_k(2,nb).ne.kv(2)).or. 462 > (brillioun_all_k(3,nb).ne.kv(3))) 463 > ierr = 1 464 end do 465 if (ierr.eq.1) then 466 if (mprint) then 467 write(*,*)"Brillioun Zone Vectors do not match!" 468 call flush(6) 469 end if 470 end if 471 else 472 if (mprint) then 473 write(*,*)"Brillioun Zone Points do not match!" 474 write(*,*)"NB = ",nbrillioun," not equal ", 475 > brillioun_nbrillioun() 476 call flush(6) 477 end if 478 ierr = 1 479 end if 480 end if 481 482 msglen = 1 483 call Parallel_Brdcst_ivalues(MASTER,msglen,ierr) 484 if (ierr.ne.0) then 485 if (taskid.eq.MASTER) call closefile(5) 486 return 487 end if 488 489 490* **** send header data to all processors **** 491 msglen = 1 492 call Parallel_Brdcst_ivalues(MASTER,msglen,version) 493 msglen = 3 494 call Parallel_Brdcst_ivalues(MASTER,msglen,nfft) 495 msglen = 9 496 call Parallel_Brdcst_values(MASTER,msglen,unita) 497 498 iatom(1) = ichar(atom(1:1)) 499 iatom(2) = ichar(atom(2:2)) 500 msglen = 2 501 call Parallel_Brdcst_ivalues(MASTER,msglen,iatom) 502 atom(1:1) = char(iatom(1)) 503 atom(2:2) = char(iatom(2)) 504 505 msglen = 1 506 call Parallel_Brdcst_ivalues(MASTER,msglen,lmax) 507 call Parallel_Brdcst_ivalues(MASTER,msglen,locp) 508 call Parallel_Brdcst_ivalues(MASTER,msglen,nbrillioun) 509 call Parallel_Brdcst_values(MASTER,msglen,rcut) 510 call Parallel_Brdcst_values(MASTER,msglen,lmbda) 511 lmmax=(lmax+1)**2 - (2*locp+1) 512 513 514 515* **** read in borb 3d blocks **** 516 borbs_tag = borb_projector_alloc(brillioun_nbrillq(),lmmax,npack1) 517 do nb=1,brillioun_nbrillioun() 518 call K1dB_ktoqp(nb,nbq,pk) 519 do n=1,lmmax 520 call C3dB_r_read(1,5,tmp2,tmp,-1,pk,.true.) 521 if (pk.eq.taskid_k) then 522 call Cram_r_pack(nbq,tmp2) 523 call borb_projector_add(borbs_tag,nbq,n,tmp2) 524 end if 525 end do 526 end do 527 528* *** close fname binary file *** 529 if (taskid.eq.MASTER) then 530 call closefile(5) 531 end if 532 533 ierr = 0 534 return 535 end 536 537* *********************************** 538* * * 539* * borbs_readall * 540* * * 541* *********************************** 542 logical function borbs_readall() 543 implicit none 544 545#include "bafdecls.fh" 546#include "borbs.fh" 547 548* **** local variables **** 549 integer ngp(3),version,nfft3d,npack1 550 integer ia,l,lm,ii,icount 551 real*8 unita(3,3) 552 integer tmp(2),tmp2(2),ierr 553 logical value,found,correct_box,value2 554 character*2 atom 555 character*4 element 556 character*50 fname 557 558* **** parallel i/o variable **** 559 integer MASTER,taskid 560 parameter(MASTER=0) 561 562* **** external functions **** 563 logical nwpw_filefind 564 integer control_ngrid 565 real*8 control_unita 566 character*12 control_boundry 567 character*4 ion_atom_qm 568 external nwpw_filefind 569 external control_ngrid 570 external control_unita 571 external control_boundry 572 external ion_atom_qm 573 integer ion_nion_qm,ion_katm_qm 574 external ion_nion_qm,ion_katm_qm 575 576 call C3dB_nfft3d(1,nfft3d) 577 call Cram_max_npack(npack1) 578 call Parallel_taskid(taskid) 579 580 value = BA_push_get(mt_dbl,(4*nfft3d),'tmp',tmp(2),tmp(1)) 581 if (.not. value) go to 9000 582 583 value = BA_push_get(mt_dbl,(2*nfft3d),'tmp2',tmp2(2),tmp2(1)) 584 if (.not. value) go to 9000 585 586* **** read pseudopotentials **** 587 do ia=1,nkatmx 588 589* **** define formatted borb name **** 590 element = ' ' 591 element = ion_atom_qm(ia) 592 l = index(element,' ') - 1 593 fname = element(1:l)//'.borb' 594 595 596 found = .false. 597 do while (.not.found) 598 599 if (nwpw_filefind(fname)) then 600 call borbs_read(fname, 601 > version, 602 > ngp,unita, 603 > atom, 604 > int_mb(lmmax(1)+ia-1), 605 > int_mb(lmax(1)+ia-1), 606 > int_mb(locp(1)+ia-1), 607 > dbl_mb(rcut(1)+ia-1), 608 > dbl_mb(lmbda(1)+ia-1), 609 > npack1, 610 > int_mb(borbs(1)+ ia-1), 611 > dbl_mb(tmp(1)),dbl_mb(tmp2(1)), 612 > ierr) 613 614 if (ierr.gt.0) then 615 value = .false. 616 go to 9000 617 end if 618 619* ************************************************************** 620* ***** logic for finding out if psp is correctly formatted **** 621* ************************************************************** 622 correct_box = .true. 623 if ( (ngp(1).ne.control_ngrid(1)) .or. 624 > (ngp(2).ne.control_ngrid(2)) .or. 625 > (ngp(3).ne.control_ngrid(3)) .or. 626 > (unita(1,1).ne.control_unita(1,1)) .or. 627 > (unita(2,1).ne.control_unita(2,1)) .or. 628 > (unita(3,1).ne.control_unita(3,1)) .or. 629 > (unita(1,2).ne.control_unita(1,2)) .or. 630 > (unita(2,2).ne.control_unita(2,2)) .or. 631 > (unita(3,2).ne.control_unita(3,2)) .or. 632 > (unita(1,3).ne.control_unita(1,3)) .or. 633 > (unita(2,3).ne.control_unita(2,3)) .or. 634 > (unita(3,3).ne.control_unita(3,3))) then 635 correct_box = .false. 636 if (taskid.eq.MASTER) then 637 write(6,*) "atomic orbitals are not correctly formatted:", 638 > fname 639 end if 640 end if 641 if (correct_box) found = .true. 642 643 end if 644 645* **** generate formatted pseudopotential atom.borb ***** 646 if (.not.found) then 647 call borbs_formatter_auto(ion_atom_qm(ia),0.0d0,0.0d0) 648 end if 649 650 end do !***do while **** 651 652 653 end do 654 655* *********************************************** 656* **** set up the index for the atomic basis **** 657* *********************************************** 658 icount = 0 659 do ii=1,ion_nion_qm() 660 ia = ion_katm_qm(ii) 661 662 do lm=1,int_mb(lmmax(1)+ia-1) 663 icount = icount + 1 664 665 int_mb(lmborb(1)+icount-1) = lm 666 int_mb(iaborb(1)+icount-1) = ia 667 int_mb(iiborb(1)+icount-1) = ii 668 int_mb(basisborb(1)+lm-1+(ii-1)*norbs_max) = icount 669 end do 670 end do 671 ibasis = icount 672 call borbs_normalize() 673 674 9000 value2 = BA_pop_stack(tmp2(2)) 675 value2 = BA_pop_stack(tmp(2)) 676 677 678 borbs_readall = value 679 return 680 end 681 682