TABLE OF CONTENTS
ABINIT/check_zarot [ Functions ]
NAME
check_zarot
FUNCTION
Debugging routine used to test zarot.
COPYRIGHT
Copyright (C) 2008-2018 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 .
INPUTS
OUTPUT
PARENTS
CHILDREN
fftbox_execute,fftbox_plan3,get_bz_item,getcprj,kdata_free,kdata_init mkkpg,paw_symcprj,pawcprj_alloc,pawcprj_free,wfd_get_cprj,wfd_sym_ur
SOURCE
27 #if defined HAVE_CONFIG_H 28 #include "config.h" 29 #endif 30 31 #include "abi_common.h" 32 33 subroutine check_zarot(npwvec,Cryst,ngfft,gvec,psps,pawang,grottb,grottbm1) 34 35 use defs_basis 36 use defs_datatypes 37 use defs_abitypes 38 use m_errors 39 use m_profiling_abi 40 41 use m_fft_mesh, only : rotate_FFT_mesh 42 use m_geometry, only : normv 43 use m_crystal, only : crystal_t 44 use m_paw_sphharm, only : initylmr 45 use m_mpinfo, only : destroy_mpi_enreg 46 use m_pawang, only : pawang_type 47 use m_pawtab, only : pawtab_type 48 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 49 50 !This section has been created automatically by the script Abilint (TD). 51 !Do not modify the following lines by hand. 52 #undef ABI_FUNC 53 #define ABI_FUNC 'check_zarot' 54 use interfaces_14_hidewrite 55 use interfaces_32_util 56 use interfaces_51_manage_mpi 57 use interfaces_56_recipspace 58 !End of the abilint section 59 60 implicit none 61 62 !Arguments ------------------------------------ 63 !scalars 64 integer,intent(in) :: npwvec 65 type(crystal_t),intent(in) :: Cryst 66 type(pawang_type),intent(in) :: pawang 67 type(pseudopotential_type),intent(in) :: psps 68 !arrays 69 integer,intent(in) :: ngfft(18) 70 integer,intent(in) :: grottb(npwvec,Cryst%timrev,Cryst%nsym),grottbm1(npwvec,Cryst%timrev,Cryst%nsym) 71 integer,intent(in) :: gvec(3,npwvec) 72 73 !Local variables------------------------------- 74 !scalars 75 integer :: aa,ig,ig_sym,iginv_sym,ii,ilpa,ilpm,isym,itim,jj,ll,lmax,mm,mqmem_ 76 integer :: nqpt_,optder,option,normchoice,npts,ix,iy,iz,ir_sym,ir !,nx,ny,nz 77 real(dp) :: err,max_diff,test,tmp,ylm_sym,rx,ry,rz 78 logical :: found !,iscompatibleFFT 79 character(len=500) :: message 80 type(MPI_type) :: Fake_MPI_enreg 81 !arrays 82 integer :: toinv(Cryst%nsym),trial(3,3),rm1(3,3) 83 integer,allocatable :: nband(:),npwarr(:),irottb(:,:) 84 real(dp),allocatable :: DS_mmpl(:,:,:),DSinv_mmpl(:,:,:),qptns(:,:),ylm_q(:,:) 85 real(dp),allocatable :: ylmgr_q(:,:,:) 86 real(dp),allocatable :: ylmr(:,:),ylmr_gr(:,:,:),nrm(:),rr(:,:) !,dum_tnons(:,:) 87 real(dp) :: search(3) 88 89 ! ************************************************************************* 90 91 write(message,'(a)')' check_zarot : enter ' 92 call wrtout(std_out,message,'COLL') 93 94 do jj=1,Cryst%nsym 95 found=.FALSE. 96 do ii=1,Cryst%nsym 97 call mati3inv(Cryst%symrec(:,:,ii),trial) 98 trial=transpose(trial) 99 if (ALL(trial==Cryst%symrec(:,:,jj))) then 100 toinv(jj)=ii 101 found=.TRUE. 102 exit 103 end if 104 end do 105 if (.not. found) then 106 MSG_ERROR("inverse not found!") 107 end if 108 end do 109 110 mqmem_=1 ; nqpt_=1 ; optder=0 111 ABI_MALLOC(npwarr,(mqmem_)) 112 ABI_MALLOC(qptns,(3,mqmem_)) 113 npwarr(:)=npwvec ; qptns(:,:)=zero 114 115 lmax=psps%mpsang-1 116 write(std_out,*)'lmax= ',lmax 117 ABI_MALLOC(ylm_q,(npwvec*mqmem_,(lmax+1)**2)) 118 ABI_MALLOC(ylmgr_q,(npwvec*mqmem_,3+6*(optder/2),(lmax+1)**2)) 119 call initmpi_seq(Fake_MPI_enreg) 120 ABI_MALLOC(nband,(1)) 121 nband=0 122 123 ! Note: dtset%nband and dtset%nsppol are not used in sequential mode 124 call initylmg(Cryst%gprimd,gvec,qptns,mqmem_,Fake_MPI_enreg,Psps%mpsang,npwvec,nband,nqpt_,npwarr,0,optder,& 125 & Cryst%rprimd,ylm_q,ylmgr_q) 126 127 call destroy_mpi_enreg(Fake_MPI_enreg) 128 129 ABI_MALLOC(DS_mmpl,(2*lmax+1,2*lmax+1,lmax+1)) 130 ABI_MALLOC(DSinv_mmpl,(2*lmax+1,2*lmax+1,lmax+1)) 131 max_diff=zero ; test=zero 132 133 do ig=1,npwvec 134 if (ig==1) cycle 135 136 do isym=1,Cryst%nsym 137 do itim=1,Cryst%timrev 138 139 ig_sym=grottb(ig,itim,isym) !index of IS G 140 DS_mmpl(:,:,:)=pawang%zarot(:,:,:,isym) 141 142 iginv_sym=grottbm1(ig,itim,isym) !index of (IS)^-1 G 143 DSinv_mmpl(:,:,:)=pawang%zarot(:,:,:,toinv(isym)) 144 145 do ll=0,lmax 146 do mm=1,2*ll+1 147 ilpm=1+ll**2+ll+(mm-1-ll) 148 ylm_sym=ylm_q(ig_sym,ilpm) !Ylm(IS G) 149 !ylm_sym=ylm_q(iginv_sym,ilpm) !Ylm(IS^-1G) 150 ! 151 ! here we calculate the symmetric 152 tmp=zero 153 do aa=1,2*ll+1 154 test=MAX(test,ABS(DS_mmpl(aa,mm,ll+1)-DSinv_mmpl(mm,aa,ll+1))) 155 ilpa=1+ll**2+ll+(aa-1-ll) 156 tmp= tmp+ ylm_q(ig,ilpa)*DS_mmpl(aa,mm,ll+1) 157 end do 158 if (itim==2) tmp=tmp*(-1)**ll 159 err=ABS(tmp-ylm_sym) !Ylm(IS G) = D_am Yma(S) (-1)**l 160 161 if (err > tol6) then 162 write(std_out,*)'WARNING check fort 77' 163 write(77,'(6(a,i3),a)')' -- ig: ',ig,' igsym: ',ig_sym,' isym ',isym,' itim:',itim,' ll: ',ll,' mm: ',(mm-1-ll)," --" 164 write(77,*)tmp,ylm_sym,ABS(tmp-ylm_sym) 165 end if 166 max_diff=MAX(max_diff,err) 167 168 end do 169 end do !itim 170 171 end do !isym 172 end do !sym 173 end do !ig 174 175 write(std_out,*)"MAX DIFF ",max_diff 176 write(std_out,*)"MAX TEST ",test 177 178 179 ABI_FREE(nband) 180 ABI_FREE(npwarr) 181 ABI_FREE(qptns) 182 ABI_FREE(ylm_q) 183 ABI_FREE(ylmgr_q) 184 185 npts = PRODUCT(ngfft(1:3)) 186 187 ABI_MALLOC(irottb,(npts,Cryst%nsym)) 188 !allocate(dum_tnons(3,Cryst%nsym)); dum_tnons=zero 189 !call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,dum_tnons,ngfft,irottb,iscompatibleFFT) 190 !if (.not.iscompatibleFFT) then 191 ! MSG_ERROR("Uncompatible FFT mesh") 192 !end if 193 !deallocate(dum_tnons) 194 195 ABI_MALLOC(rr,(3,npts)) 196 ABI_MALLOC(nrm,(npts)) 197 ii = 0 198 do iz=0,ngfft(3)-1 199 do iy=0,ngfft(2)-1 200 do ix=0,ngfft(1)-1 201 ii = ii + 1 202 if (ix <= ngfft(1)/2) then 203 rx = DBLE(ix)/ngfft(1) 204 else 205 rx = DBLE(ix-ngfft(1))/ngfft(1) 206 end if 207 if (iy <= ngfft(2)/2) then 208 ry = DBLE(iy)/ngfft(2) 209 else 210 ry = DBLE(iy-ngfft(2))/ngfft(2) 211 end if 212 if (iz <= ngfft(3)/2) then 213 rz = DBLE(iz)/ngfft(3) 214 else 215 rz = DBLE(iz-ngfft(3))/ngfft(3) 216 end if 217 rr(:,ii) = (/rx,ry,rz/) 218 nrm(ii) = normv(rr(:,ii),Cryst%rmet,"R") 219 end do 220 end do 221 end do 222 223 irottb = HUGE(0) 224 do isym=1,Cryst%nsym 225 call mati3inv(Cryst%symrel(:,:,isym),rm1) 226 rm1 = transpose(rm1) 227 do ii=1,npts 228 search = MATMUL(rm1,rr(:,ii)) 229 do jj=1,npts 230 if (ALL (ABS(search-rr(:,jj)) < tol6)) irottb(ii,isym) = jj 231 end do 232 end do 233 end do 234 235 option=1; normchoice=1 236 ABI_MALLOC(ylmr,(Psps%mpsang**2,npts)) 237 ABI_MALLOC(ylmr_gr,(3*(option/2)+6*(option/3),Psps%mpsang**2,npts)) 238 239 call initylmr(Psps%mpsang,normchoice,npts,nrm,option,rr,ylmr,ylmr_gr) 240 241 max_diff=zero ; test=zero 242 243 do isym=1,Cryst%nsym 244 do ir=1,npts 245 ir_sym = irottb(ir,isym) ! idx of R^{-1} (r-\tau) 246 if (ir_sym == HUGE(0)) then 247 write(std_out,*)"Got HUGE" 248 CYCLE 249 end if 250 251 do ll=0,lmax 252 do mm=1,2*ll+1 253 ilpm=1+ll**2+ll+(mm-1-ll) 254 ylm_sym=ylmr(ir_sym,ilpm) !Ylm(R^{-1}(r-t)) 255 !ylm_sym=ylm_q(iginv_sym,ilpm) !Ylm(IS^-1G) 256 ! 257 ! here we calculate the symmetric 258 tmp=zero 259 do aa=1,2*ll+1 260 test=MAX(test,ABS(DS_mmpl(aa,mm,ll+1)-DSinv_mmpl(mm,aa,ll+1))) 261 ilpa=1+ll**2+ll+(aa-1-ll) 262 tmp= tmp+ ylmr(ir,ilpa)*DS_mmpl(aa,mm,ll+1) 263 end do 264 !if (itim==2) tmp=tmp*(-1)**ll 265 err=ABS(tmp-ylm_sym) ! Ylm(R^{1}(r-t)) = D_am Yma(r) 266 267 if (err > tol6) then 268 write(std_out,*)'WARNING check fort 78' 269 write(77,'(5(a,i3),a)')' -- ir: ',ir,' ir_sym: ',ir_sym,' isym ',isym,' ll: ',ll,' mm: ',(mm-1-ll)," --" 270 write(77,*)tmp,ylm_sym,ABS(tmp-ylm_sym) 271 end if 272 max_diff=MAX(max_diff,err) 273 274 end do ! ll 275 end do ! mm 276 277 end do ! ir 278 end do ! isym 279 280 write(std_out,*)"MAX DIFF REAL SPACE ",max_diff 281 write(std_out,*)"MAX TEST REAL SPACE ",test 282 283 ABI_FREE(ylmr) 284 ABI_FREE(ylmr_gr) 285 ABI_FREE(irottb) 286 ABI_FREE(rr) 287 ABI_FREE(nrm) 288 289 ABI_FREE(DS_mmpl) 290 ABI_FREE(DSinv_mmpl) 291 292 end subroutine check_zarot
ABINIT/paw_check_symcprj [ Functions ]
NAME
paw_check_symcprj
FUNCTION
Test the routines used to symmetrize PAW CPRJ
COPYRIGHT
Copyright (C) 2010-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
PARENTS
sigma
CHILDREN
fftbox_execute,fftbox_plan3,get_bz_item,getcprj,kdata_free,kdata_init mkkpg,paw_symcprj,pawcprj_alloc,pawcprj_free,wfd_get_cprj,wfd_sym_ur
SOURCE
320 #if defined HAVE_CONFIG_H 321 #include "config.h" 322 #endif 323 324 #include "abi_common.h" 325 326 subroutine paw_check_symcprj(Wfd,ik_bz,band,spin,sym_mode,Cryst,Kmesh,Psps,Pawtab,Pawang,Cprj_bz) 327 328 use defs_basis 329 use defs_datatypes 330 use defs_abitypes 331 use m_errors 332 use m_profiling_abi 333 use m_fft 334 335 use m_pawang, only : pawang_type 336 use m_pawtab, only : pawtab_type 337 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 338 use m_crystal, only : crystal_t 339 use m_bz_mesh, only : kmesh_t, get_BZ_item 340 use m_wfd, only : wfd_t, wfd_get_cprj, kdata_init, kdata_free, kdata_t, wfd_sym_ur 341 342 !This section has been created automatically by the script Abilint (TD). 343 !Do not modify the following lines by hand. 344 #undef ABI_FUNC 345 #define ABI_FUNC 'paw_check_symcprj' 346 use interfaces_65_paw 347 use interfaces_66_nonlocal 348 !End of the abilint section 349 350 implicit none 351 352 !Arguments ------------------------------------ 353 !scalars 354 integer,intent(in) :: ik_bz,band,spin,sym_mode 355 type(crystal_t),intent(in) :: Cryst 356 type(kmesh_t),intent(in) :: Kmesh 357 type(Pawang_type),intent(in) :: Pawang 358 type(Pseudopotential_type),intent(in) :: Psps 359 type(wfd_t),intent(inout) :: Wfd 360 !arrays 361 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 362 type(pawcprj_type),intent(out) :: Cprj_bz(Cryst%natom,Wfd%nspinor) 363 364 !Local variables ------------------------------ 365 !scalars 366 integer :: k_sym,k_tim,ik_ibz,ig,fft_idx 367 integer :: cpopt,choice,npw_k,istwf_k,nkpg 368 integer :: iatom,iatm,isp 369 complex(dpc) :: k_eimkt 370 logical :: k_isirred 371 type(Kdata_t) :: Gdata 372 type(fftbox_plan3_t) :: plan 373 !arrays 374 integer :: k_umklp(3) 375 real(dp) :: k_bz(3) 376 real(dp),allocatable :: kpg_k(:,:),vectin(:,:) 377 complex(dpc) :: ur1_dpc(Wfd%nfft*Wfd%nspinor) 378 complex(gwpc) :: ur1(Wfd%nfft*Wfd%nspinor) 379 type(pawcprj_type),allocatable :: Cprj_srt(:,:) 380 381 !************************************************************************ 382 383 call get_BZ_item(Kmesh,ik_bz,k_bz,ik_ibz,k_sym,k_tim,k_eimkt,k_umklp,k_isirred) 384 385 if (k_isirred) then ! Symmetrization is not needed. Retrieve Cprj_ibz from Wfd and return immediately. 386 call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cprj_bz,sorted=.FALSE.) 387 RETURN 388 end if 389 390 select case (sym_mode) 391 392 case (1) ! Faster Symmetrization in reciprocal space. 393 394 call wfd_get_cprj(Wfd,band,ik_ibz,spin,Cryst,Cprj_bz,sorted=.FALSE.) 395 call paw_symcprj(ik_bz,Wfd%nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_bz) 396 397 case (2) ! Symmetrize u(r) in reciprocal space, FFT from r to G then call getcprj to obtain the symmetrized cprj. 398 399 ! Symmetrization in real space on the FFT BOX. 400 call wfd_sym_ur(Wfd,Cryst,Kmesh,band,ik_bz,spin,ur1) 401 402 istwf_k = 1 403 404 ! Init k_data associated to the G-sphere centered at k_bz. 405 call kdata_init(Gdata,Cryst,Psps,k_bz,istwf_k,Wfd%ngfft,Wfd%MPI_enreg,ecut=Wfd%ecut) 406 407 npw_k = Gdata%npw 408 ! 409 ! Compute (k+G) vectors 410 nkpg=0 411 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 412 if (nkpg>0) then 413 call mkkpg(Gdata%kg_k,kpg_k,k_bz,nkpg,npw_k) 414 end if 415 416 ABI_MALLOC(vectin,(2,npw_k*Wfd%nspinor)) 417 !ABI_CHECK(npw_k==Wfd%npwwfn,"Wrong npw") 418 ! 419 ! FFT R -> G TODO Fix issue with double precision complex. 420 ur1_dpc = ur1 421 422 call fftbox_plan3(plan,Wfd%ngfft(1:3),Wfd%ngfft(7),-1) 423 call fftbox_execute(plan,ur1_dpc) 424 425 do ig=1,npw_k ! FFT box to G-sphere. 426 fft_idx = Gdata%igfft0(ig) 427 if (fft_idx/=0) then ! G-G0 belong to the FFT mesh. 428 vectin(1,ig) = DBLE (ur1_dpc(fft_idx)) 429 vectin(2,ig) = AIMAG(ur1_dpc(fft_idx)) 430 else 431 vectin(:,ig) = zero ! Set this component to zero. 432 end if 433 end do 434 ! 435 ! Calculate SORTED cprj. 436 cpopt = 0 ! Nothing is already calculated. 437 choice = 1 438 439 ABI_DT_MALLOC(Cprj_srt,(Wfd%natom,Wfd%nspinor)) 440 call pawcprj_alloc(Cprj_srt,0,Wfd%nlmn_sort) 441 442 call getcprj(choice,cpopt,vectin,Cprj_srt,Gdata%fnl_dir0der0,& 443 & 0,Wfd%indlmn,istwf_k,Gdata%kg_k,kpg_k,k_bz,Wfd%lmnmax,Wfd%mgfft,Wfd%MPI_enreg,& 444 & Cryst%natom,Cryst%nattyp,Wfd%ngfft,Wfd%nloalg,npw_k,Wfd%nspinor,Cryst%ntypat,& 445 & Gdata%phkxred,Wfd%ph1d,Gdata%ph3d,Cryst%ucvol,1) 446 447 ABI_FREE(vectin) 448 ABI_FREE(kpg_k) 449 ! 450 ! Reorder cprj (sorted --> unsorted) 451 do iatom=1,Cryst%natom 452 iatm=Cryst%atindx(iatom) 453 do isp=1,Wfd%nspinor 454 Cprj_bz(iatom,isp)%cp=Cprj_srt(iatm,isp)%cp 455 end do 456 end do 457 call pawcprj_free(Cprj_srt) 458 ABI_DT_FREE(Cprj_srt) 459 460 call kdata_free(Gdata) 461 462 case default 463 MSG_ERROR("Wrong sym_mode") 464 end select 465 466 end subroutine paw_check_symcprj