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