1* 2* $Id$ 3* 4* $Log: not supported by cvs2svn $ 5* Revision 1.38 2009/02/08 03:26:30 bylaska 6* ...EJB 7* 8* Revision 1.37 2009/02/07 18:37:35 bylaska 9* Bassi vectorization fixes ...EJB 10* 11* Revision 1.36 2007/10/01 23:02:26 bylaska 12* removed debug io...EJB 13* 14* Revision 1.35 2007/09/29 00:34:15 bylaska 15* ...EJB 16* 17* Revision 1.34 2006/08/13 01:00:32 bylaska 18* ...EJB 19* 20* Revision 1.32 2006/01/13 02:05:26 marat 21* accelerating projector construction 22* 23* Revision 1.31 2004/10/14 21:53:29 bylaska 24* io fixes...EJB 25* 26* Revision 1.30 2004/09/04 17:56:21 bylaska 27* Added local potential to the projector file (.jpp). 28* More updates to constraint force. 29* ...EJB 30* 31* Revision 1.29 2003/10/28 19:50:51 edo 32* errquizzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz 33* 34* Revision 1.28 2003/10/21 02:05:17 marat 35* switched to new errquit by running global replace operation 36* see the script below (note it will not work on multiline errquit calls) 37* ********************************************************* 38* #!/bin/sh 39* 40* e=`find . -name "*F" -print` 41* 42* for f in $e 43* do 44* cp $f $f.bak 45* sed 's|\(^[ ].*call[ ]*errquit([^,]*\)\(,[^,]*\)\()\)|\1,0\2\3|' $f.bak > $f 46* #rm $f.bak 47* done 48* ********************************************************** 49* 50* Revision 1.27 2003/03/22 04:04:59 edo 51* removed extra arg to rayleigh_itol_... 52* 53* Revision 1.26 2003/03/05 23:16:32 bylaska 54* Commented out write statements and other minor fixes..... 55* self-consistent loop looks like it is working..... 56* ....EJB 57* 58* Revision 1.25 2003/02/13 01:58:55 bylaska 59* ...EJB 60* 61* Revision 1.24 2003/02/11 01:41:24 edo 62* eliminated f90-isms 63* 64* Revision 1.23 2003/02/09 21:27:07 marat 65* commented out some of the write statements 66* MV 67* 68* Revision 1.22 2003/02/06 06:13:09 marat 69* moved legendre and spher. harm. stuff to paw_special_functions dir 70* 71 72!************************************************** 73! 74! Name: paw_proj_init 75! 76! Purpose: initializes the paw projectors 77! 78! Created: 7/30/2002 79!************************************************** 80 subroutine paw_proj_init() 81 implicit none 82 83#include "bafdecls.fh" 84#include "paw_proj_data.fh" 85#include "paw_geom.fh" 86 87 integer paw_proj_nbasis 88 external paw_proj_nbasis 89 90* **** local variables **** 91 logical value,found 92 integer ia,npack0,npack1 93 character*4 element 94 character*50 fname 95 integer vl_ptr 96 97 98* *** get number of diff atom types *** 99 prj_nkatm = ion_nkatm() 100 101* *** number of basis functions for each atom kind *** 102 value = BA_alloc_get(mt_int,(prj_nkatm), 103 > 'prj_nbasis',prj_nbasis(2),prj_nbasis(1)) 104 value = value.and. 105 > BA_alloc_get(mt_int,(prj_nkatm), 106 > 'prj_indx',prj_indx(2),prj_indx(1)) 107 value = value.and. 108 > BA_alloc_get(mt_int,(2*prj_nkatm), 109 > 'i_prj_l',i_prj_l(2),i_prj_l(1)) 110 value = value.and. 111 > BA_alloc_get(mt_int,(2*prj_nkatm), 112 > 'i_prj_m',i_prj_m(2),i_prj_m(1)) 113 if (.not.value) call errquit('paw_proj_init: alloc heap',0,0) 114 115 116* *** generate formatted projector file if needed *** 117 call paw_proj_check_format() 118 119* *** read number of (nlm) projectors for each atom 120* from headers off the formatted files 121 call paw_proj_nbasis_read(int_mb(prj_nbasis(1))) 122 123* *** total number of projectors for diff types of atoms *** 124 prj_ntot = 0 125 do ia=1,prj_nkatm 126 prj_ntot = prj_ntot + int_mb(prj_nbasis(1)+ia-1) 127 end do 128 129 do ia=1,prj_nkatm 130 value = BA_alloc_get(mt_int,(int_mb(prj_nbasis(1)+ia-1)), 131 > 'sub_i_prj_l', 132 > int_mb(i_prj_l(1)+2*(ia-1)+1), 133 > int_mb(i_prj_l(1)+2*(ia-1))) 134 value = value.and. 135 > BA_alloc_get(mt_int,(int_mb(prj_nbasis(1)+ia-1)), 136 > 'sub_i_prj_m', 137 > int_mb(i_prj_m(1)+2*(ia-1)+1), 138 > int_mb(i_prj_m(1)+2*(ia-1))) 139 140 end do 141 142c *** allocate kspace arrays*** 143 call Pack_npack(0,npack0) 144 call Pack_npack(1,npack1) 145 value = BA_alloc_get(mt_dcpl,(prj_ntot*npack1), 146 > 'prj',prj(2),prj(1)) 147 if (.not.value) call errquit('paw_proj_init: alloc heap',0,1) 148 149* *** set prj_indx *** 150c int_mb(prj_indx(1)) = 0 151c do ia=2,prj_nkatm 152c int_mb(prj_indx(1)+ia-1) = int_mb(prj_indx(1)+ia-2) 153c > + int_mb(prj_nbasis(1) +ia-2)*npack1 154c end do 155 call paw_proj_init_sub(prj_nkatm,npack1, 156 > int_mb(prj_nbasis(1)),int_mb(prj_indx(1))) 157 158* **** allocate vloc potential space *** 159 call paw_vloc_init() 160 call paw_vloc_ptr(vl_ptr) 161 162* *** read in formatted prj's to prj common block *** 163 do ia=1,prj_nkatm 164 call paw_proj_read(ia, 165 > int_mb(int_mb(i_prj_l(1)+2*(ia-1))), 166 > int_mb(int_mb(i_prj_m(1)+2*(ia-1))), 167 > npack1, 168 > dcpl_mb(prj(1)+int_mb(prj_indx(1)+ia-1)), 169 > npack0, 170 > dbl_mb(vl_ptr+(ia-1)*npack0)) 171 > 172 end do 173 return 174 end 175 subroutine paw_proj_init_sub(n,npack,nbas,indx) 176 implicit none 177 integer n,npack,nbas(n),indx(n) 178 integer i 179 indx(1) = 0 180 do i=2,n 181 indx(i) = indx(i-1) + npack*nbas(i-1) 182 end do 183 return 184 end 185 186!************************************************** 187! 188! Name: paw_proj_i_prj 189! 190! Purpose: returns the dcpl_mb ma index of 191! the paw projectors. 192! 193! Created: 7/30/2002 194!************************************************** 195 integer function paw_proj_i_prj() 196 implicit none 197 198#include "paw_proj_data.fh" 199 200 paw_proj_i_prj = prj(1) 201 return 202 end 203 204!************************************************** 205! 206! Name: paw_proj_i_prj_atom 207! 208! Purpose: returns the dcpl_mb ma index of 209! the paw projectors. 210! 211! Created: 7/30/2002 212!************************************************** 213 integer function paw_proj_i_prj_atom(ia) 214 implicit none 215 integer ia 216 217#include "bafdecls.fh" 218#include "paw_proj_data.fh" 219 220 paw_proj_i_prj_atom = prj(1) + int_mb(prj_indx(1)+ia-1) 221 return 222 end 223 224 225!************************************************** 226! 227! Name: paw_proj_nbasis 228! 229! Purpose: returns the number of the paw projectors 230! for this kind of atom. 231! 232! Created: 7/30/2002 233!************************************************** 234 integer function paw_proj_nbasis(ia) 235 implicit none 236 integer ia 237 238#include "bafdecls.fh" 239#include "paw_proj_data.fh" 240 241 paw_proj_nbasis = int_mb(prj_nbasis(1)+ia-1) 242 return 243 end 244 245 246!************************************************** 247! 248! Name: paw_proj_total_nbasis 249! 250! Purpose: returns the number of the paw projectors 251! for this kind of atom. 252! 253! Created: 7/30/2002 254!************************************************** 255c integer function paw_proj_total_nbasis() 256c implicit none 257c integer ia 258 259c#include "bafdecls.fh" 260c#include "paw_proj_data.fh" 261c 262c paw_proj_total_nbasis = prj_ntot 263c return 264c end 265 266 267 268 269!************************************************** 270! 271! Name: paw_proj_l 272! 273! Purpose: returns the orbital quantum number 274! for the nth projector of the iath 275! kind of atom. 276! 277! Created: 7/30/2002 278!************************************************** 279 integer function paw_proj_l(n,ia) 280 implicit none 281 integer n,ia 282 283#include "bafdecls.fh" 284#include "paw_proj_data.fh" 285 286 paw_proj_l = int_mb(int_mb(i_prj_l(1)+2*(ia-1))+n-1) 287 return 288 end 289 290!************************************************** 291! 292! Name: paw_proj_m 293! 294! Purpose: returns the magnetic quantum number 295! for the nth projector of the iath 296! kind of atom. 297! 298! Created: 7/30/2002 299!************************************************** 300 integer function paw_proj_m(n,ia) 301 implicit none 302 integer n,ia 303 304#include "bafdecls.fh" 305#include "paw_proj_data.fh" 306 307 308 paw_proj_m = int_mb(int_mb(i_prj_m(1)+2*(ia-1))+n-1) 309 return 310 end 311 312 313 314!************************************************** 315! 316! Name: paw_proj_check_format 317! 318! Purpose: 319! 320! Created: 7/30/2002 321!************************************************** 322 subroutine paw_proj_check_format() 323 implicit none 324 325#include "paw_proj_data.fh" 326 327* **** local variables **** 328 logical value,found 329 integer ia,nkatm 330 character*4 element 331 character*50 fname 332 333* **** external functions **** 334 logical nwpw_filefind,paw_proj_format_ok 335 external nwpw_filefind,paw_proj_format_ok 336 337 do ia=1,prj_nkatm 338 339* **** define formatted prj name **** 340 call ion_atom_plus_suffix(ia,'.jpp',fname) 341 342* **** make sure prj names are formatted correctly **** 343 found = .false. 344 do while (.not.found) 345 if (nwpw_filefind(fname)) then 346 if (paw_proj_format_ok(fname)) found = .true. 347 end if 348 349* **** generate formatted projectors atom.jpp **** 350 if (.not.found) then 351 call paw_proj_formatter_auto(ia) 352 end if 353 end do 354 355 end do 356 357 return 358 end 359 360!************************************************** 361! Name: paw_proj_nbasis_read 362! 363! Purpose: returns the number of basis functions 364! for each of the formatted projector 365! files. 366! 367! Created: 7/30/2002 368!************************************************** 369 subroutine paw_proj_nbasis_read(nbasis) 370 implicit none 371 integer nbasis(*) 372 373#include "paw_proj_data.fh" 374 375* **** local variables **** 376 integer ia,ngrid(3) 377 real*8 unita(3,3) 378 character*50 fname 379 380 381 do ia=1,prj_nkatm 382 call ion_atom_plus_suffix(ia,'.jpp',fname) 383 call paw_proj_read_header(fname,ngrid,nbasis(ia),unita) 384 end do 385 386 return 387 end 388 389 390 391!************************************************** 392! Name: paw_proj_paw_basis_tot_nbasis 393! 394! Purpose: returns the number of basis functions 395! for each of the formatted projector 396! files. 397! 398! Created: 7/30/2002 399!************************************************** 400 integer function paw_proj_tot_nbasis() 401 implicit none 402 403#include "paw_proj_data.fh" 404 405 external paw_proj_nbasis 406 integer paw_proj_nbasis 407 408* **** local variables **** 409 integer ia 410 integer tot_nbasis 411 412* **** external functions **** 413 integer ion_natm 414 external ion_natm 415 416c !*** calculate total number of (n,l,m) projectors *** 417 tot_nbasis = 0 418c !-- loop over diff kinds of atoms -- 419 do ia=1,prj_nkatm 420c !-- ion_natm(ia) is number of atoms of kind ia -- 421 tot_nbasis = tot_nbasis 422 > + paw_proj_nbasis(ia)*ion_natm(ia) 423 end do 424 425 paw_proj_tot_nbasis = tot_nbasis 426 return 427 end 428 429!************************************************** 430! 431! Name: paw_proj_format_ok 432! 433! Purpose: returns true if header of the formatted 434! projector file agrees with control. 435! 436! Created: 7/30/2002 437!************************************************** 438 logical function paw_proj_format_ok(fname) 439 implicit none 440 character*(*) fname 441 442 443* **** local variables **** 444 logical correct_box 445 integer ngrid(3),nbasis 446 real*8 unita(3,3) 447 448* **** external functions **** 449 integer control_ngrid 450 real*8 control_unita 451 external control_ngrid 452 external control_unita 453 454 correct_box = .true. 455 call paw_proj_read_header(fname,ngrid,nbasis,unita) 456 if ( (ngrid(1).ne.control_ngrid(1)) .or. 457 > (ngrid(2).ne.control_ngrid(2)) .or. 458 > (ngrid(3).ne.control_ngrid(3)) .or. 459 > (unita(1,1).ne.control_unita(1,1)) .or. 460 > (unita(2,1).ne.control_unita(2,1)) .or. 461 > (unita(3,1).ne.control_unita(3,1)) .or. 462 > (unita(1,2).ne.control_unita(1,2)) .or. 463 > (unita(2,2).ne.control_unita(2,2)) .or. 464 > (unita(3,2).ne.control_unita(3,2)) .or. 465 > (unita(1,3).ne.control_unita(1,3)) .or. 466 > (unita(2,3).ne.control_unita(2,3)) .or. 467 > (unita(3,3).ne.control_unita(3,3))) then 468 correct_box = .false. 469 end if 470 471 paw_proj_format_ok = correct_box 472 return 473 end 474 475!************************************************** 476! 477! Name: paw_proj_read 478! 479! Purpose: read in the formatted 480! projector file. 481! 482! Created: 8/06/2002 483!************************************************** 484 subroutine paw_proj_read(ia,proj_l,proj_m,npack1,prj,npack0,vloc) 485 implicit none 486 integer ia 487 integer proj_l(*),proj_m(*) 488 integer npack1 489 complex*16 prj(npack1,*) 490 integer npack0 491 real*8 vloc(npack0) 492 493#include "bafdecls.fh" 494#include "tcgmsg.fh" 495#include "msgtypesf.h" 496 497* *** local variables *** 498 integer MASTER,taskid 499 parameter(MASTER=0) 500 501 logical value 502 integer ii,msglen,l,nbasis,ngrid(3),nfft3d 503 real*8 unita(3,3) 504 character*50 fname 505 character*255 full_filename 506 integer tmp1(2),tmp2(2) 507 complex*16 sum1 508 509 call Parallel_taskid(taskid) 510 511* **** open fname binary file **** 512 call ion_atom_plus_suffix(ia,'.jpp',fname) 513 if (taskid.eq.MASTER) then 514 515 call util_file_name_noprefix(fname,.false., 516 > .false., 517 > full_filename) 518 l = index(full_filename,' ') - 1 519 call openfile(5,full_filename,l,'r',l) 520 call iread(5,ngrid,3) 521 call dread(5,unita,9) 522 call iread(5,nbasis,1) 523 end if 524 525* **** send header data to all processors **** 526 msglen = 3 527 call BRDCST(9+MSGINT,ngrid,mitob(msglen),MASTER) 528 msglen = 9 529 call BRDCST(9+MSGDBL,unita,mdtob(msglen),MASTER) 530 msglen = 1 531 call BRDCST(9+MSGINT,nbasis,mitob(msglen),MASTER) 532 533 call D3dB_nfft3d(1,nfft3d) 534 value = BA_push_get(mt_dcpl,nfft3d,'tmp1',tmp1(2),tmp1(1)) 535 value = value.and. 536 > BA_push_get(mt_dcpl,nfft3d,'tmp2',tmp2(2),tmp2(1)) 537 if (.not.value) call errquit('paw_proj_read: push stack',0,0) 538 539* **** read in projectors **** 540 do ii=1,nbasis 541 if (taskid.eq.MASTER) then 542 call iread(5,proj_l(ii),1) 543 call iread(5,proj_m(ii),1) 544 end if 545 call BRDCST(9+MSGINT,proj_l(ii),mitob(msglen),MASTER) 546 call BRDCST(9+MSGINT,proj_m(ii),mitob(msglen),MASTER) 547 call D3dB_c_read(1,5,dcpl_mb(tmp1(1)),dcpl_mb(tmp2(1)),-1) 548 call Pack_c_pack(1,dcpl_mb(tmp1(1))) 549 call Pack_c_Copy(1,dcpl_mb(tmp1(1)),prj(1,ii)) 550 end do 551 552* **** read in local potential **** 553 call D3dB_t_read(1,5,dcpl_mb(tmp1(1)),dcpl_mb(tmp2(1)),-1) 554 call Pack_t_pack(0,dcpl_mb(tmp1(1))) 555 call Pack_t_Copy(0,dcpl_mb(tmp1(1)),vloc) 556 557 value = BA_pop_stack(tmp2(2)) 558 value = value.and.BA_pop_stack(tmp1(2)) 559 if (.not.value) call errquit('paw_proj_read: pop stack',0,1) 560 561 if (taskid.eq.MASTER) then 562 call closefile(5) 563 end if 564 565 return 566 end 567 568 569 570!************************************************** 571! 572! Name: paw_proj_read_header 573! 574! Purpose: read in the header of the formatted 575! projector file. 576! 577! Created: 7/30/2002 578!************************************************** 579 subroutine paw_proj_read_header(fname,ngrid,nbasis,unita) 580 implicit none 581 character*(*) fname 582 integer ngrid(3),nbasis 583 real*8 unita(3,3) 584 585#include "tcgmsg.fh" 586#include "msgtypesf.h" 587 588* *** local variables *** 589 integer MASTER,taskid 590 parameter(MASTER=0) 591 integer msglen,l 592 character*255 full_filename 593 594 call Parallel_taskid(taskid) 595 596* **** open fname binary file **** 597 if (taskid.eq.MASTER) then 598 call util_file_name_noprefix(fname,.false., 599 > .false., 600 > full_filename) 601 l = index(full_filename,' ') - 1 602 call openfile(5,full_filename,l,'r',l) 603 call iread(5,ngrid,3) 604 call dread(5,unita,9) 605 call iread(5,nbasis,1) 606 call closefile(5) 607 end if 608 609* **** send header data to all processors **** 610 msglen = 3 611 call BRDCST(9+MSGINT,ngrid,mitob(msglen),MASTER) 612 msglen = 9 613 call BRDCST(9+MSGDBL,unita,mdtob(msglen),MASTER) 614 msglen = 1 615 call BRDCST(9+MSGINT,nbasis,mitob(msglen),MASTER) 616 617 return 618 end 619 620 621 622!************************************************** 623! 624! Name: paw_proj_formatter_auto 625! 626! Purpose: read in the header of the formatted 627! projector file. 628! 629! Created: 7/30/2002 630!************************************************** 631 subroutine paw_proj_formatter_auto(ia) 632 implicit none 633 integer ia 634 635#include "bafdecls.fh" 636#include "bessel_transform.fh" 637#include "paw_params.fh" 638#include "paw_basis.fh" 639 640 641 !*** local variables *** 642 integer MASTER,taskid 643 parameter (MASTER=0) 644 645 real*8 small 646 parameter (small=1.0d-9) 647 648 logical value 649 integer ii,i,npack1,nfft3d,nr,basis_nbasis,l,m,nbasis,ngrid(3),jj 650 integer i_rgrid,i_prj,ps_ptr 651 real*8 unita(3,3),gg,log_amesh 652 integer tmp(2),prj(2),rayleigh(2),Gx(2),Gy(2),Gz(2),f(2) 653 integer g(2),gm(2),ng,ray(2) 654 655 character*50 jppname,atomname 656 character*255 full_filename 657 658 !*** external functions *** 659 integer G_indx 660 real*8 lattice_unita,lattice_omega 661 external G_indx 662 external lattice_unita,lattice_omega 663 664 665 call ion_atom_plus_suffix(ia,'_basis',atomname) 666 call ion_atom_plus_suffix(ia,'.jpp',jppname) 667 668 669 !*** read in projectors from _basis file *** 670 i_rgrid = paw_basis_i_rgrid(ia) 671 i_prj = paw_basis_i_prj_ps(ia) 672 ps_ptr = paw_basis_i_v_ps(ia) 673 basis_nbasis = int_mb(paw_basis_i_nbasis(ia)) 674 nr = int_mb(paw_basis_i_ngrid(ia)) 675 log_amesh = dbl_mb(paw_basis_i_log_amesh(ia)) 676 677 678 nbasis = 0 679 do ii=1,basis_nbasis 680 l = int_mb(paw_basis_i_orb_l(ia)+ii-1) 681 nbasis = nbasis + 2*l+1 682 end do 683 684 685 call Parallel_taskid(taskid) 686 unita(1,1) = lattice_unita(1,1) 687 unita(2,1) = lattice_unita(2,1) 688 unita(3,1) = lattice_unita(3,1) 689 unita(1,2) = lattice_unita(1,2) 690 unita(2,2) = lattice_unita(2,2) 691 unita(3,2) = lattice_unita(3,2) 692 unita(1,3) = lattice_unita(1,3) 693 unita(2,3) = lattice_unita(2,3) 694 unita(3,3) = lattice_unita(3,3) 695 call D3dB_nx(1,ngrid(1)) 696 call D3dB_ny(1,ngrid(2)) 697 call D3dB_nz(1,ngrid(3)) 698 699 700* **** open jppname binary file and write header **** 701 if (taskid.eq.MASTER) then 702 call util_file_name_noprefix(jppname,.false., 703 > .false., 704 > full_filename) 705 l = index(full_filename,' ') - 1 706 call openfile(6,full_filename,l,'w',l) 707 call iwrite(6,ngrid,3) 708 call dwrite(6,unita,9) 709 call iwrite(6,nbasis,1) 710 end if 711 712* **** compute bessel transforms **** 713 call D3dB_nfft3d(1,nfft3d) 714 call Pack_npack(1,npack1) 715 value = BA_push_get(mt_dbl,nfft3d, 716 > 'rayleigh',rayleigh(2),rayleigh(1)) 717 value = value.and. 718 > BA_push_get(mt_dcpl,nfft3d, 719 > 'prj',prj(2),prj(1)) 720 value = value.and. 721 > BA_push_get(mt_dcpl,nfft3d, 722 > 'tmp',tmp(2),tmp(1)) 723 value = value.and. 724 > BA_push_get(mt_dbl,nr,'f',f(2),f(1)) 725 value = value.and. 726 > BA_push_get(mt_int,nfft3d,'gm',gm(2),gm(1)) 727 value = value.and. 728 > BA_push_get(mt_dbl,nfft3d,'ray small',ray(2),ray(1)) 729 value = value.and. 730 > BA_push_get(mt_dbl,nfft3d,'g',g(2),g(1)) 731 if (.not.value) 732 > call errquit('paw_proj_formatter_auto, push stack',0,0) 733 734 Gx(1) = G_indx(1) 735 Gy(1) = G_indx(2) 736 Gz(1) = G_indx(3) 737 738 call dfill(nfft3d,0.0d0,dbl_mb(ray(1)),1) 739 call dfill(nfft3d,-1.0d0,dbl_mb(g(1)),1) 740 call ifill(nfft3d,0,int_mb(gm(1)),1) 741c 742 ng = 0 743 do i=1,nfft3d 744 gg = dsqrt( 745 > dbl_mb(Gx(1)+i-1)*dbl_mb(Gx(1)+i-1) 746 > + dbl_mb(Gy(1)+i-1)*dbl_mb(Gy(1)+i-1) 747 > + dbl_mb(Gz(1)+i-1)*dbl_mb(Gz(1)+i-1) ) 748 do ii=1,ng 749 if(abs(gg-dbl_mb(g(1)+ii-1)).lt.1.0d-8) then 750 int_mb(gm(1)+i-1) = ii 751 goto 1 752 end if 753 end do 754 ng = ng + 1 755 dbl_mb(g(1)+ng-1) = gg 756 int_mb(gm(1)+i-1) = ng 7571 continue 758 end do 759 760* ***** format projectors ***** 761 do ii=1,basis_nbasis 762 l = int_mb(paw_basis_i_orb_l(ia)+ii-1) 763 764 765 do i=1,ng 766 767 gg=dbl_mb(g(1)+i-1) 768 if (gg.gt.small) then 769 dbl_mb(ray(1)+i-1) 770 > = spher_bessel_transform(l,nr,log_amesh, 771 > dbl_mb(i_rgrid), 772 > dbl_mb(i_prj+(ii-1)*nr), 773 > gg) 774 else 775 dbl_mb(ray(1)+i-1) 776 > = spher_bessel0_transform(l,nr,log_amesh, 777 > dbl_mb(i_rgrid), 778 > dbl_mb(i_prj+(ii-1)*nr)) 779 end if 780 end do 781 782c do i=1,nfft3d 783c dbl_mb(rayleigh(1)+i-1) = 784c > dbl_mb(ray(1)+int_mb(gm(1)+i-1)-1) 785c end do 786 call paw_proj_sub2(nfft3d,int_mb(gm(1)), 787 > dbl_mb(ray(1)),dbl_mb(rayleigh(1))) 788 789 do m=-l,l 790* *** generate Ylm *** 791 call spher_harmonics_generate(l,m,nfft3d, 792 > dbl_mb(Gx(1)), 793 > dbl_mb(Gy(1)), 794 > dbl_mb(Gz(1)), 795 > dcpl_mb(prj(1))) 796 797* *** multiply Ylm and spherical besse 798c call D3dB_tc_Mul(1,dbl_mb(rayleigh(1)), 799c > dcpl_mb(prj(1)), 800c > dcpl_mb(prj(1))) 801 call D3dB_tc_Mul2(1,dbl_mb(rayleigh(1)), 802 > dcpl_mb(prj(1))) 803 804 call rayleigh_itol_scaling(l,nfft3d, 805 > dcpl_mb(prj(1))) 806 807 if (taskid.eq.MASTER) then 808 call iwrite(6,l,1) 809 call iwrite(6,m,1) 810 end if 811 !call Pack_c_unpack(1,dcpl_mb(prj(1))) 812 call D3dB_c_write(1,6,dcpl_mb(prj(1)),dcpl_mb(tmp(1)),0) 813 end do ! m loop 814 end do ! basis_nbasis loop 815 816 817 818* ***** format local potential ***** 819c do i=1,nr 820c dbl_mb(f(1)+i-1) = dbl_mb(i_rgrid+i-1) 821c > *dbl_mb(ps_ptr+i-1) 822c end do 823 call paw_proj_sub1(nr,dbl_mb(i_rgrid),dbl_mb(ps_ptr), 824 > dbl_mb(f(1))) 825 826 do i=1,ng 827 gg=dbl_mb(g(1)+i-1) 828 829 if (gg.gt.small) then 830 dbl_mb(ray(1)+i-1) 831 > = fourpi 832 > *spher_bessel_transform(0,nr,log_amesh, 833 > dbl_mb(i_rgrid), 834 > dbl_mb(f(1)), 835 > gg) 836 else 837 dbl_mb(ray(1)+i-1) 838 > = fourpi 839 > *spher_bessel0_transform(0,nr,log_amesh, 840 > dbl_mb(i_rgrid), 841 > dbl_mb(f(1))) 842 end if 843 844 end do 845c do i=1,nfft3d 846c dbl_mb(rayleigh(1)+i-1) = 847c > dbl_mb(ray(1)+int_mb(gm(1)+i-1)-1) 848c end do 849 call paw_proj_sub2(nfft3d,int_mb(gm(1)), 850 > dbl_mb(ray(1)),dbl_mb(rayleigh(1))) 851 852 call D3dB_t_write(1,6,dbl_mb(rayleigh(1)),dcpl_mb(tmp(1)),0) 853 854 855 if (taskid.eq.MASTER) then 856 call closefile(6) 857 write(6,*) 858 l = index(full_filename,' ') - 1 859 write(6,*) " Generated formatted filename:",full_filename(1:l) 860 end if 861 862 value = BA_pop_stack(g(2)) 863 value = value.and.BA_pop_stack(ray(2)) 864 value = value.and.BA_pop_stack(gm(2)) 865 value = value.and.BA_pop_stack(f(2)) 866 value = value.and.BA_pop_stack(tmp(2)) 867 value = value.and.BA_pop_stack(prj(2)) 868 value = value.and.BA_pop_stack(rayleigh(2)) 869 if (.not.value) 870 > call errquit('paw_proj_formatter_auto, pop stack',1,0) 871 872 return 873 end 874 875 subroutine paw_proj_sub1(n,a,b,c) 876 implicit none 877 integer n 878 real*8 a(n),b(n),c(n) 879 integer i 880 do i=1,n 881 c(i) = a(i)*b(i) 882 end do 883 return 884 end 885 886 subroutine paw_proj_sub2(n,indx,ray,rayl) 887 implicit none 888 integer n,indx(n) 889 real*8 ray(*),rayl(*) 890 integer i 891 do i=1,n 892 rayl(i) = ray(indx(i)) 893 end do 894 return 895 end 896 897!************************************************** 898! 899! Name: rayleigh_itol_scaling 900! 901! Purpose: This routine should be rewritten 902! to improve performance. 903! This routine scales Y <-- (-i)**l * Y 904! 905! Created: 8/06/2002 906!************************************************** 907 subroutine rayleigh_itol_scaling(l,nfft3d,Y) 908 implicit none 909 integer l,nfft3d 910 complex*16 Y(*) 911 912 real*8 lattice_omega 913 external lattice_omega 914#include "paw_params.fh" 915 916 !*** local variables *** 917 integer k 918 complex*16 coef 919 920 if (mod(l,4).eq.0) then 921 coef = dcmplx(1.0d0,0.0d0) 922 else if (mod(l,4).eq.1) then 923 coef = dcmplx(0.0d0,-1.0d0) 924 else if (mod(l,4).eq.2) then 925 coef = dcmplx(-1.0d0,0.0d0) 926 else if (mod(l,4).eq.3) then 927 coef = dcmplx(0.0d0,1.0d0) 928 end if 929 930 coef = coef*fourpi/sqrt(lattice_omega()) 931 do k=1,nfft3d 932 Y(k) = coef*Y(k) 933 end do 934 return 935 end 936 937 subroutine paw_proj_weghts_set() 938 939 return 940 end 941!************************************************** 942! 943! Name: paw_proj_end 944! 945! Purpose: removes space used by the paw projectors 946! 947! Created: 7/30/2002 948!************************************************** 949 subroutine paw_proj_end() 950 implicit none 951 952#include "bafdecls.fh" 953#include "paw_proj_data.fh" 954 955 !*** local variables *** 956 logical value 957 integer ia 958 959 value = .true. 960 do ia=1,prj_nkatm 961 value = value.and.BA_free_heap(int_mb(i_prj_l(1)+2*(ia-1)+1)) 962 value = value.and.BA_free_heap(int_mb(i_prj_m(1)+2*(ia-1)+1)) 963 end do 964 value = value.and.BA_free_heap(i_prj_l(2)) 965 value = value.and.BA_free_heap(i_prj_m(2)) 966 value = value.and.BA_free_heap(prj_indx(2)) 967 value = value.and.BA_free_heap(prj_nbasis(2)) 968 value = value.and.BA_free_heap(prj(2)) 969 if (.not.value) call errquit('paw_proj_end: dealloc heap',0,0) 970 971 972 !*** deallocate local potential *** 973 call paw_vloc_end() 974 975 return 976 end 977 978 979 980