TABLE OF CONTENTS
- ABINIT/cgtk_change_gsphere
- ABINIT/cgtk_fixphase
- ABINIT/cgtk_rotate
- ABINIT/cgtk_rotate_symrec
- ABINIT/m_cgtk
ABINIT/cgtk_change_gsphere [ Functions ]
NAME
cgtk_change_gsphere
FUNCTION
Transfer the G components of ndat wavefunctions from one sphere to another one. Can also be used to change the value of istwfk e.g. 2 --> 1
INPUTS
ndat = Number of wavefunctions to transform. npw1, npw2 = Number of plane-waves in the (input, output) G-sphere istwf1, istwf2 = Storage mode of (input, output) wavefunctions. kg1(3,npw1), kg2(3,npw2) = Input/Output G-sphere cg1(2,npw1,ndat) = Input wavefunctions on kg1 sphere with istwf1 mode. work_ngfft(18)=Specify work dimensions. Must be large enough to accomodate kg1 and kg2
OUTPUT
cg2(2,npw2,ndat) = Output wavefunctions on kg2 sphere with istwf2 mode. work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) = Workspace array
SOURCE
408 subroutine cgtk_change_gsphere(ndat, npw1, istwf1, kg1, cg1, npw2, istwf2, kg2, cg2, work_ngfft, work) 409 410 !Arguments ------------------------------------ 411 !scalars 412 integer,intent(in) :: ndat,npw1,npw2,istwf1,istwf2 413 !arrays 414 integer,intent(in) :: kg1(3,npw1),kg2(3,npw2) 415 integer,intent(in) :: work_ngfft(18) 416 real(dp),intent(inout) :: cg1(2,npw1,ndat) ! TODO: Should be intent(in) but need to change sphere 417 real(dp),intent(out) :: cg2(2,npw2,ndat) 418 real(dp),intent(out) :: work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) 419 420 !Local variables ------------------------------ 421 !scalars 422 integer :: n1,n2,n3,n4,n5,n6,idat 423 424 !************************************************************************ 425 426 n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3) 427 n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6) 428 429 !print *, "npw1", npw1, "npw2", npw2 430 431 do idat=1,ndat 432 ! Insert cg1 in work array taking into account istwf1 (intent in) 433 call sphere(cg1(:,:,idat),1,npw1,work,n1,n2,n3,n4,n5,n6,kg1,istwf1,to_box,me_g0,no_shift,identity_3d,one) 434 435 ! Extract cg2 from work array taking into account istwf2 436 call sphere(cg2(:,:,idat),1,npw2,work,n1,n2,n3,n4,n5,n6,kg2,istwf2,to_sph,me_g0,no_shift,identity_3d,one) 437 end do 438 439 end subroutine cgtk_change_gsphere
ABINIT/cgtk_fixphase [ Functions ]
NAME
cgtk_fixphase
FUNCTION
Fix phase of all bands. Keep normalization but maximize real part (minimize imag part). Also fix the sign of real part by setting the first non-zero element to be positive. See also fxphas_seq in m_cgtools
INPUTS
cg(2,mcg)= contains the wavefunction |c> coefficients. gsc(2,mgsc)= if useoverlap==1, contains the S|c> coefficients, where S is an overlap matrix. icg=shift to be applied on the location of data in the array cg igsc=shift to be applied on the location of data in the array gsc istwfk=input option parameter that describes the storage of wfs (set to 1 if usual complex vectors) mcg=size of second dimension of cg mgsc=size of second dimension of gsc mpi_enreg=information about MPI parallelization nband_k=number of bands npw_k=number of planewaves useoverlap=describe the overlap of wavefunctions: 0: no overlap (S=Identi0,ty_matrix) 1: wavefunctions are overlapping
OUTPUT
cg(2,mcg)=same array with altered phase. gsc(2,mgsc)= same array with altered phase.
SOURCE
473 subroutine cgtk_fixphase(cg, gsc, icg, igsc, istwfk, mcg, mgsc, mpi_enreg, nband_k, npw_k, useoverlap, cprj, nspinor) 474 475 !Arguments ------------------------------------ 476 !scalars 477 integer,intent(in) :: icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap 478 type(MPI_type),intent(in) :: mpi_enreg 479 !arrays 480 real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc*useoverlap) 481 type(pawcprj_type),intent(inout),optional,target :: cprj(:,:) 482 integer,intent(in),optional :: nspinor 483 484 !Local variables------------------------------- 485 !scalars 486 logical :: do_cprj 487 integer :: iband,ierr,ii,indx,ncprj 488 real(dp) :: cim,cre,gscim,gscre,quotient,root1,root2,saa,sab,sbb,theta,thppi,xx,yy 489 character(len=500) :: msg 490 !arrays 491 real(dp) :: buffer2(nband_k,2),buffer3(nband_k,3),tsec(2) 492 real(dp),allocatable :: cimb(:),creb(:),saab(:),sabb(:),sbbb(:) !,sarr(:,:) 493 494 ! ************************************************************************* 495 496 do_cprj=.false. 497 if (present(cprj)) then 498 do_cprj=.true. 499 ncprj = size(cprj,2) 500 if (ncprj/=nband_k*nspinor) then 501 ABI_ERROR('bad size for cprj') 502 end if 503 end if 504 505 !The general case, where a complex phase indeterminacy is present 506 if(istwfk==1)then 507 508 ABI_MALLOC(cimb,(nband_k)) 509 ABI_MALLOC(creb,(nband_k)) 510 ABI_MALLOC(saab,(nband_k)) 511 ABI_MALLOC(sabb,(nband_k)) 512 ABI_MALLOC(sbbb,(nband_k)) 513 cimb(:)=zero ; creb(:)=zero 514 515 ! Loop over bands 516 ! TODO: MG store saa arrays in sarr(3,nband_k) to reduce false sharing. 517 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nband_k,icg,npw_k,cg,saab,sbbb,sabb) 518 do iband=1,nband_k 519 indx=icg+(iband-1)*npw_k 520 521 ! Compute several sums over Re, Im parts of c 522 saa=zero; sbb=zero; sab=zero 523 do ii=1+indx,npw_k+indx 524 saa=saa+cg(1,ii)*cg(1,ii) 525 sbb=sbb+cg(2,ii)*cg(2,ii) 526 sab=sab+cg(1,ii)*cg(2,ii) 527 end do 528 saab(iband)=saa 529 sbbb(iband)=sbb 530 sabb(iband)=sab 531 end do 532 533 ! XG030513 : MPIWF : should transmit saab,sbbb,sabb from the procs 534 ! of the WF group to the master processor of the WF group 535 if (mpi_enreg%paral_kgb == 1) then 536 buffer3(:,1)=saab(:) 537 buffer3(:,2)=sbbb(:) 538 buffer3(:,3)=sabb(:) 539 call timab(48,1,tsec) 540 call xmpi_sum(buffer3,mpi_enreg%comm_fft,ierr) 541 if (mpi_enreg%paral_spinor==1) then 542 call xmpi_sum(buffer3,mpi_enreg%comm_spinor,ierr) 543 end if 544 call timab(48,2,tsec) 545 saab(:)=buffer3(:,1) 546 sbbb(:)=buffer3(:,2) 547 sabb(:)=buffer3(:,3) 548 end if 549 550 ! XG030513 : MPIWF this loop should only be executed by the master of the WF group 551 552 if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then 553 do iband=1,nband_k 554 indx=icg+(iband-1)*npw_k 555 556 saa=saab(iband) 557 sbb=sbbb(iband) 558 sab=sabb(iband) 559 560 ! Get phase angle theta 561 if (sbb+saa>tol8)then 562 if(abs(sbb-saa)>tol8*(sbb+saa) .or. 2*abs(sab)>tol8*(sbb+saa))then 563 if (abs(sbb-saa)>tol8*abs(sab)) then 564 quotient=sab/(sbb-saa) 565 theta=0.5_dp*atan(2.0_dp*quotient) 566 else 567 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 568 theta=0.25_dp*(pi-(sbb-saa)/sab) 569 end if 570 ! Check roots to get theta for max Re part 571 root1=cos(theta)**2*saa+sin(theta)**2*sbb-2.0_dp*cos(theta)*sin(theta)*sab 572 thppi=theta+0.5_dp*pi 573 root2=cos(thppi)**2*saa+sin(thppi)**2*sbb-2.0_dp*cos(thppi)*sin(thppi)*sab 574 if (root2>root1) theta=thppi 575 else 576 ! The real part vector and the imaginary part vector are orthogonal, and of same norm. Strong indeterminacy. 577 ! Will determine the first non-zero coefficient, and fix its phase 578 ! Hypothesis : there is at least one non-zero element on the master node ... 579 do ii=1+indx,npw_k+indx 580 cre=cg(1,ii) 581 cim=cg(2,ii) 582 if(cre**2+cim**2>tol8**2*(saa+sbb))then 583 if(cre**2>tol8**2**cim**2)then 584 theta=atan(cim/cre) 585 else 586 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 587 theta=pi/2-cre/cim 588 end if 589 exit 590 end if 591 end do 592 end if 593 else 594 write(msg,'(a,i0,5a)')& 595 & 'The eigenvector with band ',iband,' has zero norm.',ch10,& 596 & 'This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)',ch10,& 597 & 'Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut' 598 ABI_ERROR(msg) 599 end if 600 601 xx=cos(theta) 602 yy=sin(theta) 603 604 ! Here, set the first non-zero element to be positive 605 ! Comment the next nine lines to recover the behaviour of pre v3.1.3g 606 ! Hypothesis : there is at least one non-zero element on the master node ... 607 do ii=1+indx,npw_k+indx 608 cre=cg(1,ii) 609 cim=cg(2,ii) 610 cre=xx*cre-yy*cim 611 if(abs(cre)>tol8)exit 612 end do 613 if(cre<zero)then 614 xx=-xx ; yy=-yy 615 end if 616 617 creb(iband)=xx 618 cimb(iband)=yy 619 620 end do 621 end if 622 623 ! XG030513 : MPIWF : should transmit creb(:),cimb(:) of the master 624 ! processor of the WF group to the others procs of the WF group 625 if (mpi_enreg%paral_kgb == 1) then 626 call timab(48,1,tsec) 627 buffer2(:,1)=creb(:) 628 buffer2(:,2)=cimb(:) 629 call xmpi_sum(buffer2,mpi_enreg%comm_fft,ierr) 630 if (mpi_enreg%paral_spinor==1) then 631 call xmpi_sum(buffer2,mpi_enreg%comm_spinor,ierr) 632 end if 633 call timab(48,2,tsec) 634 creb(:)=buffer2(:,1) 635 cimb(:)=buffer2(:,2) 636 end if 637 638 ! MG TODO: Scaling can be done with zscal 639 !$OMP PARALLEL DO PRIVATE(indx,xx,yy,cre,cim,gscre,gscim) 640 do iband=1,nband_k 641 indx=icg+(iband-1)*npw_k 642 643 xx=creb(iband) 644 yy=cimb(iband) 645 ! Alter phase of array |cg> 646 do ii=1+indx,npw_k+indx 647 cre=cg(1,ii) 648 cim=cg(2,ii) 649 cg(1,ii)=xx*cre-yy*cim 650 cg(2,ii)=xx*cim+yy*cre 651 end do 652 if (do_cprj) call pawcprj_zaxpby((/zero,zero/),(/xx,yy/),cprj(:,nspinor*(iband-1)+1:nspinor*iband),& 653 & cprj(:,nspinor*(iband-1)+1:nspinor*iband)) 654 655 ! Alter phase of array S|cg> 656 if (useoverlap==1) then 657 indx=igsc+(iband-1)*npw_k 658 do ii=1+indx,npw_k+indx 659 gscre=gsc(1,ii) 660 gscim=gsc(2,ii) 661 gsc(1,ii)=xx*gscre-yy*gscim 662 gsc(2,ii)=xx*gscim+yy*gscre 663 end do 664 end if 665 end do ! iband 666 667 ABI_FREE(cimb) 668 ABI_FREE(creb) 669 ABI_FREE(saab) 670 ABI_FREE(sabb) 671 ABI_FREE(sbbb) 672 673 else ! if istwfk/=1. Storages that take into account the time-reversal symmetry : the freedom is only a sign freedom 674 675 ABI_MALLOC(creb,(nband_k)) 676 creb(:)=zero 677 ! XG030513 : MPIWF : this loop should be done only by the master processor of the WF group 678 679 if (mpi_enreg%paral_kgb==0.or.mpi_enreg%me_fft==0) then 680 681 ! Loop over bands 682 do iband=1,nband_k 683 684 indx=icg+(iband-1)*npw_k 685 686 ! Here, set the first non-zero real element to be positive 687 do ii=1+indx,npw_k+indx 688 cre=cg(1,ii) 689 if(abs(cre)>tol8)exit 690 end do 691 creb(iband)=cre 692 693 end do ! iband 694 695 end if 696 ! XG030513 : MPIWF : should transmit cre(:) of the master processor of the WF group to the others 697 if (mpi_enreg%paral_kgb == 1) then 698 call timab(48,1,tsec) 699 call xmpi_sum(creb,mpi_enreg%comm_fft,ierr) 700 if (mpi_enreg%paral_spinor==1) then 701 call xmpi_sum(creb,mpi_enreg%comm_spinor,ierr) 702 end if 703 call timab(48,2,tsec) 704 end if 705 706 do iband=1,nband_k 707 cre=creb(iband) 708 if(cre<zero)then 709 indx=icg+(iband-1)*npw_k 710 do ii=1+indx,npw_k+indx 711 cg(1,ii)=-cg(1,ii) 712 cg(2,ii)=-cg(2,ii) 713 end do 714 if (do_cprj) call pawcprj_zaxpby((/zero,zero/),(/-one,zero/),cprj(:,iband:iband),cprj(:,iband:iband)) 715 if(useoverlap==1)then 716 indx=igsc+(iband-1)*npw_k 717 do ii=1+indx,npw_k+indx 718 gsc(1,ii)=-gsc(1,ii) 719 gsc(2,ii)=-gsc(2,ii) 720 end do 721 end if 722 end if 723 end do ! iband 724 725 ABI_FREE(creb) 726 end if ! istwfk 727 728 end subroutine cgtk_fixphase
ABINIT/cgtk_rotate [ Functions ]
NAME
cgtk_rotate
FUNCTION
Reconstruct wavefunction cg2 in the BZ from the symmetrical image cg1 by applying a symmetry operation. Note that there are two possible conventions for mapping k-points: 1) k2 = T symrel(:,:, isym)^t k1 + g0 (note transpose of symrel) 2) k2 = T symrec(:,:, isym) k1 + g0 where T is for time-reversal (itimrev) This routine assumes the FIRST convention that, unfortunately, is not very handy. The second convention, indeed, is the most natural one when mapping k-points.
INPUTS
cryst=crystalline structure kpt1(3)=k-point in cg1. isym=Index of symmetry operation (symrel^T convention) itimrev=1 if time-reversal is needed else 0. g0(3)=g0 vector nspinor=Number of spinor components. ndat=Number of wavefunctions npw1, npw2=Number of G-vectors in kg1 and kg2. kg1(3,npw1), kg2(3,npw2) = G vectors in cg1, and cg2. istwf1, istwf2= Storage mode for cg1 and cg2 work_ngfft(18)= Specifies the size of the workspace array work. IMPORTANT: must be large enough to accoung for all possible shifts of the g-sphere. The caller is responsible for computing the max size needed to handle all the possible symmetrization. cg1(2, npw1, nspinor, ndat)=Wavefunctions in the IBZ
OUTPUT
cg2(2, npw2, nspinor, ndat)= symmetrized wavefunctions. work(2, work_ngfft(4), work_ngfft(5), work_ngfft(6))) = workspace array. See comments in INPUTS section.
NOTES
Inspired to wfconv.
SOURCE
104 subroutine cgtk_rotate(cryst, kpt1, isym, itimrev, g0, nspinor, ndat, & 105 npw1, kg1, npw2, kg2, istwf1, istwf2, cg1, cg2, work_ngfft, work) 106 107 !Arguments ------------------------------------ 108 !scalars 109 integer,intent(in) :: isym, itimrev, nspinor, ndat, npw1, npw2, istwf1, istwf2 110 type(crystal_t),intent(in) :: cryst 111 !arrays 112 integer,intent(in) :: g0(3), kg1(3,npw1), kg2(3,npw2), work_ngfft(18) 113 real(dp),intent(in) :: kpt1(3), cg1(2,npw1,nspinor,ndat) 114 real(dp),intent(out) :: cg2(2,npw2,nspinor,ndat) 115 real(dp),intent(out) :: work(2,work_ngfft(4),work_ngfft(5),work_ngfft(6)) !*ndat) for threads? 116 117 !Local variables ------------------------------ 118 !scalars 119 integer :: n1,n2,n3,n4,n5,n6,ipw,idat,isp 120 real(dp) :: arg,ar,ai,bi,br,spinrots,spinrotx,spinroty,spinrotz 121 logical :: have_phase 122 !arrays 123 integer,parameter :: atindx(1) = 1 124 integer :: symrec(3,3), symrel(3,3) 125 real(dp) :: phktnons(2,1), tau(3), spinrot(4), tsec(2) 126 real(dp),allocatable :: phase1d(:,:), phase3d(:,:), wavef1(:,:) 127 128 !************************************************************************ 129 130 ! Keep track of total time spent. 131 call timab(1780, 1, tsec) 132 133 ABI_CHECK_IRANGE(itimrev, 0, 1, "itimrev should be in [0, 1]") 134 135 n1 = work_ngfft(1); n2 = work_ngfft(2); n3 = work_ngfft(3) 136 n4 = work_ngfft(4); n5 = work_ngfft(5); n6 = work_ngfft(6) 137 138 symrel = cryst%symrel(:,:,isym) 139 call mati3inv(symrel, symrec) ! symrec = symrel^{-1t} 140 tau = cryst%tnons(:,isym) 141 have_phase = sum(tau ** 2) > tol8 142 143 ! Compute rotation in spinor space 144 if (nspinor == 2) call getspinrot(cryst%rprimd, spinrot, symrel) 145 if (itimrev == 1) symrec = -symrec 146 147 ! Need to compute phase factors associated with nonsymmorphic translations? 148 if (have_phase) then 149 150 ! Although the routine getph is originally written for atomic phase factors, it does precisely what we want 151 ABI_MALLOC(phase1d, (2, (2*n1+1)+(2*n2+1)+(2*n3+1))) 152 call getph(atindx, 1, n1, n2, n3, phase1d, tau) 153 154 arg = two_pi * (kpt1(1)*tau(1) + kpt1(2)*tau(2) + kpt1(3)*tau(3)) 155 phktnons(1, 1) = cos(arg) 156 phktnons(2, 1) = sin(arg) 157 ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k1 + G).tnons ) 158 ABI_MALLOC(phase3d, (2, npw1)) 159 call ph1d3d(1, 1, kg1, 1, 1, npw1, n1, n2, n3, phktnons, phase1d, phase3d) 160 ABI_FREE(phase1d) 161 end if 162 163 ABI_MALLOC(wavef1, (2, npw1)) 164 165 do idat=1,ndat 166 do isp=1,nspinor 167 wavef1 = cg1(:,:,isp,idat) 168 169 if (have_phase) then 170 ! Multiply by phase factors due to nonsymmorphic translations. 171 do ipw=1,npw1 172 ar = phase3d(1,ipw) * wavef1(1,ipw) - phase3d(2,ipw) * wavef1(2,ipw) 173 ai = phase3d(2,ipw) * wavef1(1,ipw) + phase3d(1,ipw) * wavef1(2,ipw) 174 wavef1(1, ipw) = ar 175 wavef1(2, ipw) = ai 176 end do 177 end if 178 179 ! Take into account time-reversal symmetry for SCALAR wavefunctions, if needed. 180 if (itimrev == 1 .and. nspinor == 1) wavef1(2, :npw1) = -wavef1(2, :npw1) 181 182 ! Insert wavef1 in work array. 183 call sphere(wavef1,ndat1,npw1,work,n1,n2,n3,n4,n5,n6,kg1,istwf1,to_box,me_g0,no_shift,identity_3d,one) 184 185 ! Apply rotation + g0 and extract data on the kg2 sphere: cg2(g) = work(S(g + g0)) 186 call sphere(cg2(:,:,isp,idat),ndat1,npw2,work,n1,n2,n3,n4,n5,n6,kg2,istwf2,to_sph,me_g0,g0,symrec,one) 187 end do ! isp 188 189 if (nspinor == 2) then 190 if (itimrev == 1) then 191 ! Take care of time-reversal symmetry, if needed 192 ! 1) Exchange spin-up and spin-down. 193 ! 2) Make complex conjugate of one component, and change sign of other component 194 do ipw=1,npw2 195 ! Here, change sign of real part 196 ar = -cg2(1,ipw,1,idat) 197 ai = cg2(2,ipw,1,idat) 198 ! Here, change sign of imaginary part 199 cg2(1,ipw,1,idat) = cg2(1,ipw,2,idat) 200 cg2(2,ipw,1,idat) = -cg2(2,ipw,2,idat) 201 cg2(1,ipw,2,idat) = ar 202 cg2(2,ipw,2,idat) = ai 203 end do 204 end if ! itimrev==1 205 206 ! Rotation in spinor space (see also wfconv) 207 spinrots = spinrot(1); spinrotx = spinrot(2); spinroty = spinrot(3); spinrotz = spinrot(4) 208 do ipw=1,npw2 209 ar = cg2(1,ipw,1,idat) 210 ai = cg2(2,ipw,1,idat) 211 br = cg2(1,ipw,2,idat) 212 bi = cg2(2,ipw,2,idat) 213 cg2(1,ipw,1,idat) = spinrots*ar - spinrotz*ai + spinroty*br - spinrotx*bi 214 cg2(2,ipw,1,idat) = spinrots*ai + spinrotz*ar + spinroty*bi + spinrotx*br 215 cg2(1,ipw,2,idat) = -spinroty*ar - spinrotx*ai + spinrots*br + spinrotz*bi 216 cg2(2,ipw,2,idat) = -spinroty*ai + spinrotx*ar + spinrots*bi - spinrotz*br 217 end do 218 end if 219 end do ! idat 220 221 ABI_FREE(wavef1) 222 ABI_SFREE(phase3d) 223 224 call timab(1780, 2, tsec) 225 226 end subroutine cgtk_rotate
ABINIT/cgtk_rotate_symrec [ Functions ]
NAME
cgtk_rotate_symrec
FUNCTION
Reconstruct wavefunction cg2 in the BZ from the symmetrical image cg1 by applying a symmetry operation. Note that there are two possible conventions for mapping k-points: 1) k2 = T symrel(:,:, isym)^t k1 + g0 (note transpose of symrel) 2) k2 = T symrec(:,:, isym) k1 + g0 where T is for time-reversal (itimrev) This routine assumes the SECOND convention. For scalar wavefunctions, we have (with S being a symrec operation) 1) u_{Sk}(g) = e^{-i(Sk + g).tau)} u_k(S^{-1} g) if g0 = 0 and no TR 2) u_{-k}(g) = u_{k}(-g)^* for TR 3) u_{k+g0}(g) = u_{k}(g+g0) if g0 != 0
ABINIT/m_cgtk [ Modules ]
NAME
m_cgtk
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_cgtk 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_xmpi 28 use m_time 29 30 use m_fstrings, only : itoa, sjoin 31 use defs_abitypes, only : MPI_type 32 use m_symtk, only : mati3inv 33 use m_geometry, only : getspinrot 34 use m_crystal, only : crystal_t 35 use m_fftcore, only : sphere 36 use m_kg, only : ph1d3d, getph 37 use m_pawcprj, only : pawcprj_type, pawcprj_zaxpby 38 39 implicit none 40 41 private