TABLE OF CONTENTS
ABINIT/nonlop_test [ Functions ]
NAME
nonlop_test
FUNCTION
This routine is supposed to be used only for testing purpose. It tests the "nonlop" routine application (non-local operator) with respect to Finite Differences. It is not supposed to be used standalone, but via the nonlop_dfpt_test.py script to be found in ~abinit/scripts/post_processing/nonlop_dfpt_test directory. This Python script launches Abinit (several datasets) and analyse the result, in order to compare <Psi_i|H^(i)|Psi_j> compute with DFPT or Finite Differences. H^(i) is the ith derivative of the Hamiltonian with respect to one or several perturbation(s).
COPYRIGHT
Copyright (C) 2017-2018 ABINIT group (MT) 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
cg(2,mcg)=wavefunctions (may be read from disk file) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) istwfk(nkpt)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis kpt(3,nkpt)=k points in reduced coordinates mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mkmem=number of k points treated by this node. mpi_enreg=informations about MPI parallelization mpw= maximum number of plane waves my_natom=number of atoms treated by current processor natom=number of atoms in cell. nband(nkpt)=number of bands at each k point nfft=number of FFT grid points ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt=number of k points in Brillouin zone nloalg(3)=governs the choice of the algorithm for non-local operator. npwarr(nkpt)=number of planewaves in basis and boundary at each k nspden=Number of spin Density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(natom)=type of each atom xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
PARENTS
afterscfloop
CHILDREN
destroy_hamiltonian,dotprod_g,init_hamiltonian,initylmg load_k_hamiltonian,load_spin_hamiltonian,mkffnl,mkkpg,nonlop
SOURCE
65 #if defined HAVE_CONFIG_H 66 #include "config.h" 67 #endif 68 69 #include "abi_common.h" 70 71 subroutine nonlop_test(cg,eigen,istwfk,kg,kpt,mband,mcg,mgfft,mkmem,mpi_enreg,mpw,my_natom,natom,& 72 & nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,ntypat,& 73 & paw_ij,pawtab,ph1d,psps,rprimd,typat,xred) 74 75 use defs_basis 76 use defs_datatypes 77 use defs_abitypes 78 use m_profiling_abi 79 use m_xmpi 80 use m_errors 81 use m_hamiltonian 82 use m_pawtab 83 use m_paw_ij 84 use m_pawcprj 85 use m_cgtools 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'nonlop_test' 91 use interfaces_32_util 92 use interfaces_56_recipspace 93 use interfaces_66_nonlocal, except_this_one => nonlop_test 94 !End of the abilint section 95 96 implicit none 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer :: mband,mcg,mgfft,mkmem,mpw,my_natom,natom,nfft,nkpt,nspden,nspinor,nsppol,ntypat 101 type(MPI_type),intent(inout) :: mpi_enreg 102 type(pseudopotential_type),intent(in) :: psps 103 !arrays 104 integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol),ngfft(18),nloalg(3),npwarr(nkpt),typat(natom) 105 real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol),kpt(3,nkpt), ph1d(2,3*(2*mgfft+1)*natom) 106 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 107 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 108 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 109 !Local variables------------------------------- 110 !scalars 111 integer,parameter :: ndtset_test=6,tim_nonlop=4 112 integer,save :: idtset=0 113 integer :: bandpp,bdtot_index,blocksize,choice,cplex,cpopt,dimffnl,iatom,iatom_only 114 integer :: iband,iband_last,iband_test,iblock,icg,ider_ffnl,idir,idir_ffnl,idir_nonlop 115 integer :: ii,ikg,ikpt,ilm,inlout,isppol,istwf_k,me_distrb,my_nspinor,nband_k,nblockbd 116 integer :: nkpg,nnlout,npw_k,paw_opt,signs,spaceComm 117 logical :: ex 118 character(len=100) :: strg 119 real(dp) :: argr,argi 120 type(gs_hamiltonian_type) :: gs_hamk 121 !arrays 122 integer,allocatable :: kg_k(:,:) 123 real(dp) :: kpoint(3),rmet(3,3) 124 real(dp),allocatable :: cwavef(:,:),cwavef_out(:,:),enl(:,:,:),enlout(:),kpg_k(:,:),lambda(:) 125 real(dp),allocatable :: scwavef_out(:,:),ylm(:,:),ylmgr(:,:,:),ylm_k(:,:),ylmgr_k(:,:,:) 126 real(dp),allocatable,target :: ffnl(:,:,:,:),ph3d(:,:,:) 127 type(pawcprj_type) :: cwaveprj(1,1) 128 129 !************************************************************************* 130 131 !Increment dataset counter 132 idtset=idtset+1 133 if (idtset<=2) return 134 135 !Data from parallelism 136 spaceComm=mpi_enreg%comm_kpt 137 me_distrb=mpi_enreg%me_kpt 138 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 139 140 !Initialize Hamiltonian datastructure 141 call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,& 142 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,& 143 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 144 & mpi_spintab=mpi_enreg%my_isppoltab,paw_ij=paw_ij,ph1d=ph1d) 145 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 146 147 !Check for existence of files in the current directory\ 148 ! and set parameters accordingly 149 choice=1 ; idir=0 ; signs=1 150 if(idtset<ndtset_test)then 151 inquire(file='config/signs1',exist=ex) ; if(ex) signs=1 152 inquire(file='config/signs2',exist=ex) ; if(ex) signs=2 153 do ii=1,100 154 if (ii< 10) write(unit=strg,fmt='(a13,i1)') "config/choice",ii 155 if (ii>=10) write(unit=strg,fmt='(a13,i2)') "config/choice",ii 156 inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii 157 end do 158 do ii=1,9 159 write(unit=strg,fmt='(a11,i1)') "config/idir",ii 160 inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii 161 end do 162 else 163 inquire(file='config/signsdfpt1',exist=ex) ; if(ex)signs=1 164 inquire(file='config/signsdfpt2',exist=ex) ; if(ex)signs=2 165 do ii=1,100 166 if (ii< 10) write(unit=strg,fmt='(a17,i1)') "config/choicedfpt",ii 167 if (ii>=10) write(unit=strg,fmt='(a17,i2)') "config/choicedfpt",ii 168 inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii 169 end do 170 do ii=1,36 171 if (ii< 10) write(unit=strg,fmt='(a15,i1)') "config/idirdfpt",ii 172 if (ii>=10) write(unit=strg,fmt='(a15,i2)') "config/idirdfpt",ii 173 inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii 174 end do 175 end if 176 iatom=1 ; iband_test=-1 177 do ii=1,50 178 if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iatom",ii 179 if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iatom",ii 180 inquire(file=trim(strg),exist=ex) ; if(ex) iatom=ii 181 if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iband",ii 182 if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iband",ii 183 inquire(file=trim(strg),exist=ex) ; if(ex) iband_test=ii 184 end do 185 186 !Set parameters for the "nonlop" routine according to users choice 187 cpopt=-1 ; paw_opt=3*psps%usepaw 188 inquire(file='config/dij',exist=ex);if(ex) paw_opt=1*psps%usepaw 189 if(signs==1)then 190 iatom_only=-1 ; idir_ffnl=0 191 idir_nonlop=0 ; cplex=1 192 if(choice==1)then 193 ider_ffnl=0 194 nnlout=1 ; inlout=1 195 end if 196 if(choice==2)then 197 ider_ffnl=0 198 nnlout=3*natom ; inlout=3*(iatom-1)+idir 199 end if 200 if(choice==3)then 201 ider_ffnl=1 202 nnlout=6 ; inlout=idir 203 end if 204 if(choice==5)then 205 ider_ffnl=1 206 nnlout=3 ; inlout=idir 207 end if 208 if(choice==51.or.choice==52)then 209 ider_ffnl=1 ; cplex=2 210 nnlout=6 ; inlout=2*idir-1 211 end if 212 if(choice==54)then 213 ider_ffnl=2 ; cplex=2 214 nnlout=18*natom ; inlout=18*(iatom-1)+2*idir-1 215 end if 216 if(choice==55)then 217 ider_ffnl=2 ; cplex=2 218 nnlout=36 ; inlout=2*idir-1 219 end if 220 if(choice==8)then 221 ider_ffnl=2 222 nnlout=6 ; inlout=idir 223 end if 224 if(choice==81)then 225 ider_ffnl=2 ; cplex=2 226 nnlout=18 ; inlout=2*idir-1 227 end if 228 else if(signs==2)then 229 nnlout=1 ; inlout =1 ; cplex=1 230 idir_nonlop=idir ; iatom_only=-1 231 if(choice==1)then 232 ider_ffnl=0 ; idir_ffnl=0 233 end if 234 if(choice==2)then 235 iatom_only=iatom 236 ider_ffnl=0 ; idir_ffnl=0 237 end if 238 if(choice==3)then 239 ider_ffnl=1 ; idir_ffnl=-7 240 end if 241 if(choice==5)then 242 ider_ffnl=1 ; idir_ffnl=4 243 end if 244 if(choice==51.or.choice==52)then 245 ider_ffnl=1 ; idir_ffnl=4 ; cplex=2 246 end if 247 if(choice==54)then 248 iatom_only=iatom 249 ider_ffnl=2 ; idir_ffnl=4 ; cplex=2 250 end if 251 if(choice==8)then 252 ider_ffnl=2 ; idir_ffnl=4 253 end if 254 if(choice==81)then 255 ider_ffnl=2 ; idir_ffnl=4 ; cplex=2 256 end if 257 end if 258 259 !Set parameters for the "mkffnl" routine according to users choice 260 dimffnl=1+ider_ffnl 261 if (ider_ffnl==1.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=2+2*psps%useylm 262 if (ider_ffnl==2.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=3+7*psps%useylm 263 if (ider_ffnl==1.and.idir_ffnl==-7) dimffnl=2+5*psps%useylm 264 if (idir_ffnl>-7.and.idir_ffnl<0) dimffnl=2 265 266 !Write recognizable statement in log file 267 write(std_out,'(2(a,i2),(a,i1),2(a,i2),(a,i1),(a,i2),(a,i1),2(a,i2))') & 268 & "TESTDFPT: choice=",choice,", idir(mkffnl)=",idir_ffnl,& 269 & ", ider(mkffnl)=",ider_ffnl,", dimffnl=",dimffnl,& 270 & ", idir(nonlop)=",idir_nonlop,", signs=",signs,& 271 & ", iatom=",iatom_only,", paw_opt=",paw_opt,& 272 & ", nnlout=",nnlout,", inlout=",inlout 273 274 !Compute all spherical harmonics and gradients 275 ABI_ALLOCATE(ylm,(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 276 ABI_ALLOCATE(ylmgr,(mpw*mkmem,9,psps%mpsang*psps%mpsang*psps%useylm)) 277 if (psps%useylm==1) then 278 call initylmg(gs_hamk%gprimd,kg,kpt,mkmem,mpi_enreg,psps%mpsang,mpw,nband,nkpt,& 279 & npwarr,nsppol,2,rprimd,ylm,ylmgr) 280 else 281 ylm=zero ; ylmgr=zero 282 end if 283 284 !No loop over spins; only do the first one 285 bdtot_index=0 ; icg=0 ; isppol=1 286 287 !Continue to initialize the Hamiltonian (PAW DIJ coefficients) 288 call load_spin_hamiltonian(gs_hamk,isppol,with_nonlocal=.true.) 289 290 !No loop over k points; only do the first one 291 ikg=0 ; ikpt=1 292 293 nband_k=nband(ikpt+(isppol-1)*nkpt) 294 istwf_k=istwfk(ikpt) 295 npw_k=npwarr(ikpt) 296 kpoint(:)=kpt(:,ikpt) 297 298 !My spin/kpoint or not? 299 if(.not.proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 300 301 ! Parallelism over FFT and/or bands: define sizes and tabs 302 bandpp=mpi_enreg%bandpp 303 nblockbd=nband_k/bandpp 304 blocksize=nband_k/nblockbd 305 306 ! Several allocations 307 ABI_ALLOCATE(lambda,(blocksize)) 308 ABI_ALLOCATE(enlout,(nnlout*blocksize)) 309 ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize)) 310 ABI_ALLOCATE(cwavef_out,(2,npw_k)) 311 if (paw_opt>=3) then 312 ABI_ALLOCATE(scwavef_out,(2,npw_k)) 313 ABI_ALLOCATE(enl,(0,0,0)) 314 else 315 ABI_ALLOCATE(scwavef_out,(0,0)) 316 ABI_ALLOCATE(enl,(gs_hamk%dimekb1,gs_hamk%dimekb2,gs_hamk%nspinor**2)) 317 enl(:,:,:)=one 318 end if 319 320 ! Compute (k+G) vectors and associated spherical harmonics 321 nkpg=3*nloalg(3) 322 ABI_ALLOCATE(kg_k,(3,mpw)) 323 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 324 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 325 ABI_ALLOCATE(ylmgr_k,(npw_k,9,psps%mpsang*psps%mpsang*psps%useylm)) 326 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 327 if (nkpg>0) then 328 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 329 end if 330 if (psps%useylm==1) then 331 do ilm=1,psps%mpsang*psps%mpsang 332 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 333 end do 334 if (ider_ffnl>=1) then 335 do ilm=1,psps%mpsang*psps%mpsang 336 do ii=1,3+6*(ider_ffnl/2) 337 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 338 end do 339 end do 340 end if 341 end if 342 343 ! Compute non-local form factors 344 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 345 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,& 346 & gs_hamk%gmet,gs_hamk%gprimd,ider_ffnl,idir_ffnl,psps%indlmn,kg_k,kpg_k,& 347 & gs_hamk%kpt_k,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 348 & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,& 349 & psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 350 351 ! Load k-dependent part in the Hamiltonian datastructure 352 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk)) 353 call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,& 354 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.) 355 356 do iblock=1,nblockbd 357 358 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k) 359 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle 360 361 ! Select a specific band or all 362 if (iband==iband_test.or.iband_test==-1) then 363 364 ! Load contribution from block(n,k) 365 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 366 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 367 lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 368 369 ! Call NONLOP 370 if (paw_opt<3) then 371 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,& 372 & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,& 373 & iatom_only=iatom_only,enl=enl) 374 else 375 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,& 376 & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,& 377 & iatom_only=iatom_only) 378 end if 379 380 ! Post-processing if nonlop is called with specific options 381 if (signs==2) then 382 if (paw_opt<3) then 383 call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,cwavef_out,& 384 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 385 else 386 call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,scwavef_out,& 387 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 388 end if 389 enlout(inlout)=argr 390 end if 391 if (signs==1.and.choice==1) then 392 call dotprod_g(argr,argi,istwf_k,npw_k,1,cwavef,cwavef,& 393 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 394 enlout(:)=enlout(:)+argr 395 end if 396 397 ! Write recognizable statements in log file 398 if (idtset<ndtset_test) then 399 write(std_out,'(a,i3,es24.16)') "TESTDFPT_df: ",iband,enlout(inlout) 400 else 401 write(std_out,'(a,i3,es24.16)') "TESTDFPT_dfpt:",iband,enlout(inlout) 402 end if 403 404 end if 405 406 end do ! End of loop on block of bands 407 408 ! Increment indexes (not used here because only one spin/kpoint) 409 icg=icg+npw_k*my_nspinor*nband_k 410 ikg=ikg+npw_k 411 412 end if ! Not my spin/kpoint 413 414 bdtot_index=bdtot_index+nband_k 415 416 !Memory deallocations 417 ABI_DEALLOCATE(enl) 418 ABI_DEALLOCATE(enlout) 419 ABI_DEALLOCATE(lambda) 420 ABI_DEALLOCATE(ph3d) 421 ABI_DEALLOCATE(ffnl) 422 ABI_DEALLOCATE(cwavef) 423 ABI_DEALLOCATE(cwavef_out) 424 ABI_DEALLOCATE(scwavef_out) 425 ABI_DEALLOCATE(kg_k) 426 ABI_DEALLOCATE(kpg_k) 427 ABI_DEALLOCATE(ylm_k) 428 ABI_DEALLOCATE(ylmgr_k) 429 ABI_DEALLOCATE(ylm) 430 ABI_DEALLOCATE(ylmgr) 431 call destroy_hamiltonian(gs_hamk) 432 433 end subroutine nonlop_test