TABLE OF CONTENTS
ABINIT/m_nonlop_test [ Modules ]
NAME
m_nonlop_test
FUNCTION
COPYRIGHT
Copyright (C) 2017-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_nonlop_test 22 23 implicit none 24 25 private
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).
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=information 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
SOURCE
83 subroutine nonlop_test(cg,eigen,istwfk,kg,kpt,mband,mcg,mgfft,mkmem,mpi_enreg,mpw,my_natom,natom,& 84 & nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,ntypat,& 85 & paw_ij,pawtab,ph1d,psps,rprimd,typat,xred) 86 87 use defs_basis 88 use m_abicore 89 use m_xmpi 90 use m_errors 91 use m_hamiltonian 92 use m_pawtab 93 use m_paw_ij 94 use m_pawcprj 95 use m_cgtools 96 97 98 use defs_datatypes, only : pseudopotential_type 99 use defs_abitypes, only : MPI_type 100 use m_kg, only : mkkpg 101 use m_initylmg, only : initylmg 102 use m_mkffnl, only : mkffnl 103 use m_mpinfo, only : proc_distrb_cycle 104 use m_nonlop, only : nonlop 105 106 !Arguments ------------------------------------ 107 !scalars 108 integer :: mband,mcg,mgfft,mkmem,mpw,my_natom,natom,nfft,nkpt,nspden,nspinor,nsppol,ntypat 109 type(MPI_type),intent(inout) :: mpi_enreg 110 type(pseudopotential_type),intent(in) :: psps 111 !arrays 112 integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol),ngfft(18),nloalg(3),npwarr(nkpt),typat(natom) 113 real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol),kpt(3,nkpt), ph1d(2,3*(2*mgfft+1)*natom) 114 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 115 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 116 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 117 !Local variables------------------------------- 118 !scalars 119 integer,parameter :: ndtset_test=6,tim_nonlop=4 120 integer,save :: idtset=0 121 integer :: bandpp,bdtot_index,blocksize,choice,cplex,cpopt,dimffnl,iatm,iatom,iatom_only 122 integer :: iband,iband_last,iband_test,iblock,icg,ider_ffnl,idir,idir_ffnl,idir_nonlop 123 integer :: ii,ikg,ikpt,ilm,inlout,isppol,istwf_k,me_distrb,my_nspinor,nband_k,nblockbd 124 integer :: nkpg,nnlout,npw_k,paw_opt,signs,spaceComm 125 logical :: ex 126 character(len=100) :: strg 127 real(dp) :: argr,argi 128 type(gs_hamiltonian_type) :: gs_hamk 129 !arrays 130 integer,allocatable :: kg_k(:,:) 131 real(dp) :: kpoint(3),rmet(3,3) 132 real(dp),allocatable :: cwavef(:,:),cwavef_out(:,:),enl(:,:,:,:),enlout(:),kpg_k(:,:),lambda(:) 133 real(dp),allocatable :: scwavef_out(:,:),ylm(:,:),ylmgr(:,:,:),ylm_k(:,:),ylmgr_k(:,:,:) 134 real(dp),allocatable,target :: ffnl(:,:,:,:),ph3d(:,:,:) 135 type(pawcprj_type) :: cwaveprj(1,1) 136 137 !************************************************************************* 138 139 !Increment dataset counter 140 idtset=idtset+1 141 if (idtset<=2) return 142 143 !Data from parallelism 144 spaceComm=mpi_enreg%comm_kpt 145 me_distrb=mpi_enreg%me_kpt 146 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 147 148 !Initialize Hamiltonian datastructure 149 call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,& 150 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,& 151 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 152 & mpi_spintab=mpi_enreg%my_isppoltab,paw_ij=paw_ij,ph1d=ph1d) 153 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 154 155 !Check for existence of files in the current directory\ 156 ! and set parameters accordingly 157 choice=1 ; idir=0 ; signs=1 158 if(idtset<ndtset_test)then 159 inquire(file='config/signs1',exist=ex) ; if(ex) signs=1 160 inquire(file='config/signs2',exist=ex) ; if(ex) signs=2 161 do ii=1,100 162 if (ii< 10) write(unit=strg,fmt='(a13,i1)') "config/choice",ii 163 if (ii>=10) write(unit=strg,fmt='(a13,i2)') "config/choice",ii 164 inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii 165 end do 166 do ii=1,9 167 write(unit=strg,fmt='(a11,i1)') "config/idir",ii 168 inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii 169 end do 170 else 171 inquire(file='config/signsdfpt1',exist=ex) ; if(ex)signs=1 172 inquire(file='config/signsdfpt2',exist=ex) ; if(ex)signs=2 173 do ii=1,100 174 if (ii< 10) write(unit=strg,fmt='(a17,i1)') "config/choicedfpt",ii 175 if (ii>=10) write(unit=strg,fmt='(a17,i2)') "config/choicedfpt",ii 176 inquire(file=trim(strg),exist=ex) ; if(ex) choice=ii 177 end do 178 do ii=1,36 179 if (ii< 10) write(unit=strg,fmt='(a15,i1)') "config/idirdfpt",ii 180 if (ii>=10) write(unit=strg,fmt='(a15,i2)') "config/idirdfpt",ii 181 inquire(file=trim(strg),exist=ex) ; if(ex) idir=ii 182 end do 183 end if 184 iatom=1 ; iband_test=-1 185 do ii=1,50 186 if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iatom",ii 187 if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iatom",ii 188 inquire(file=trim(strg),exist=ex) ; if(ex) iatom=ii 189 if (ii< 10) write(unit=strg,fmt='(a12,i1)') "config/iband",ii 190 if (ii>=10) write(unit=strg,fmt='(a12,i2)') "config/iband",ii 191 inquire(file=trim(strg),exist=ex) ; if(ex) iband_test=ii 192 end do 193 194 !Set parameters for the "nonlop" routine according to users choice 195 cpopt=-1 ; paw_opt=3*psps%usepaw ; iatm=gs_hamk%atindx(iatom) 196 inquire(file='config/dij',exist=ex);if(ex) paw_opt=1*psps%usepaw 197 if(signs==1)then 198 iatom_only=-1 ; idir_ffnl=0 199 idir_nonlop=0 ; cplex=1 200 if(choice==1)then 201 ider_ffnl=0 202 nnlout=1 ; inlout=1 203 end if 204 if(choice==2)then 205 ider_ffnl=0 206 nnlout=3*natom ; inlout=3*(iatm-1)+idir ! Atoms are type-sorted in enlout() 207 end if 208 if(choice==3)then 209 ider_ffnl=1 210 nnlout=6 ; inlout=idir 211 end if 212 if(choice==5)then 213 ider_ffnl=1 214 nnlout=3 ; inlout=idir 215 end if 216 if(choice==51.or.choice==52)then 217 ider_ffnl=1 ; cplex=2 218 nnlout=6 ; inlout=2*idir-1 219 end if 220 if(choice==54)then 221 ider_ffnl=2 ; cplex=2 222 nnlout=18*natom ; inlout=18*(iatm-1)+2*idir-1 ! Atoms are type-sorted in enlout() 223 end if 224 if(choice==55)then 225 ider_ffnl=2 ; cplex=2 226 nnlout=36 ; inlout=2*idir-1 227 end if 228 if(choice==8)then 229 ider_ffnl=2 230 nnlout=6 ; inlout=idir 231 end if 232 if(choice==81)then 233 ider_ffnl=2 ; cplex=2 234 nnlout=18 ; inlout=2*idir-1 235 end if 236 else if(signs==2)then 237 nnlout=1 ; inlout =1 ; cplex=1 238 idir_nonlop=idir ; iatom_only=-1 239 if(choice==1)then 240 ider_ffnl=0 ; idir_ffnl=0 241 end if 242 if(choice==2)then 243 iatom_only=iatom 244 ider_ffnl=0 ; idir_ffnl=0 245 end if 246 if(choice==3)then 247 ider_ffnl=1 ; idir_ffnl=-7 248 end if 249 if(choice==5)then 250 ider_ffnl=1 ; idir_ffnl=4 251 end if 252 if(choice==51.or.choice==52)then 253 ider_ffnl=1 ; idir_ffnl=4 ; cplex=2 254 end if 255 if(choice==54)then 256 iatom_only=iatom 257 ider_ffnl=2 ; idir_ffnl=4 ; cplex=2 258 end if 259 if(choice==8)then 260 ider_ffnl=2 ; idir_ffnl=4 261 end if 262 if(choice==81)then 263 ider_ffnl=2 ; idir_ffnl=4 ; cplex=2 264 end if 265 end if 266 267 !Set parameters for the "mkffnl" routine according to users choice 268 dimffnl=1+ider_ffnl 269 if (ider_ffnl==1.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=2+2*psps%useylm 270 if (ider_ffnl==2.and.(idir_ffnl==0.or.idir_ffnl==4)) dimffnl=3+7*psps%useylm 271 if (ider_ffnl==1.and.idir_ffnl==-7) dimffnl=2+5*psps%useylm 272 if (idir_ffnl>-7.and.idir_ffnl<0) dimffnl=2 273 274 !Write recognizable statement in log file 275 write(std_out,'(2(a,i2),(a,i1),2(a,i2),(a,i1),(a,i2),(a,i1),2(a,i2))') & 276 & "TESTDFPT: choice=",choice,", idir(mkffnl)=",idir_ffnl,& 277 & ", ider(mkffnl)=",ider_ffnl,", dimffnl=",dimffnl,& 278 & ", idir(nonlop)=",idir_nonlop,", signs=",signs,& 279 & ", iatom=",iatom_only,", paw_opt=",paw_opt,& 280 & ", nnlout=",nnlout,", inlout=",inlout 281 282 !Compute all spherical harmonics and gradients 283 ABI_MALLOC(ylm,(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 284 ABI_MALLOC(ylmgr,(mpw*mkmem,9,psps%mpsang*psps%mpsang*psps%useylm)) 285 if (psps%useylm==1) then 286 call initylmg(gs_hamk%gprimd,kg,kpt,mkmem,mpi_enreg,psps%mpsang,mpw,nband,nkpt,& 287 & npwarr,nsppol,2,rprimd,ylm,ylmgr) 288 else 289 ylm=zero ; ylmgr=zero 290 end if 291 292 !No loop over spins; only do the first one 293 bdtot_index=0 ; icg=0 ; isppol=1 294 295 !Continue to initialize the Hamiltonian (PAW DIJ coefficients) 296 call gs_hamk%load_spin(isppol,with_nonlocal=.true.) 297 298 !No loop over k points; only do the first one 299 ikg=0 ; ikpt=1 300 301 nband_k=nband(ikpt+(isppol-1)*nkpt) 302 istwf_k=istwfk(ikpt) 303 npw_k=npwarr(ikpt) 304 kpoint(:)=kpt(:,ikpt) 305 306 !My spin/kpoint or not? 307 if(.not.proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 308 309 ! Parallelism over FFT and/or bands: define sizes and tabs 310 bandpp=mpi_enreg%bandpp 311 nblockbd=nband_k/bandpp 312 blocksize=nband_k/nblockbd 313 314 ! Several allocations 315 ABI_MALLOC(lambda,(blocksize)) 316 ABI_MALLOC(enlout,(nnlout*blocksize)) 317 ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize)) 318 ABI_MALLOC(cwavef_out,(2,npw_k)) 319 if (paw_opt>=3) then 320 ABI_MALLOC(scwavef_out,(2,npw_k)) 321 ABI_MALLOC(enl,(0,0,0,0)) 322 else 323 ABI_MALLOC(scwavef_out,(0,0)) 324 ABI_MALLOC(enl,(gs_hamk%dimekb1,gs_hamk%dimekb2,gs_hamk%nspinor**2,1)) 325 enl(:,:,:,:)=one 326 end if 327 328 ! Compute (k+G) vectors and associated spherical harmonics 329 nkpg=3*nloalg(3) 330 ABI_MALLOC(kg_k,(3,mpw)) 331 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 332 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 333 ABI_MALLOC(ylmgr_k,(npw_k,9,psps%mpsang*psps%mpsang*psps%useylm)) 334 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 335 if (nkpg>0) then 336 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 337 end if 338 if (psps%useylm==1) then 339 do ilm=1,psps%mpsang*psps%mpsang 340 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 341 end do 342 if (ider_ffnl>=1) then 343 do ilm=1,psps%mpsang*psps%mpsang 344 do ii=1,3+6*(ider_ffnl/2) 345 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 346 end do 347 end do 348 end if 349 end if 350 351 ! Compute non-local form factors 352 ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 353 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,& 354 & gs_hamk%gmet,gs_hamk%gprimd,ider_ffnl,idir_ffnl,psps%indlmn,kg_k,kpg_k,& 355 & gs_hamk%kpt_k,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 356 & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,& 357 & psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 358 359 ! Load k-dependent part in the Hamiltonian datastructure 360 ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk)) 361 call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,& 362 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.) 363 364 do iblock=1,nblockbd 365 366 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k) 367 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle 368 369 ! Select a specific band or all 370 if (iband==iband_test.or.iband_test==-1) then 371 372 ! Load contribution from block(n,k) 373 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 374 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 375 lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 376 377 ! Call NONLOP 378 if (paw_opt<3) then 379 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,& 380 & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,& 381 & iatom_only=iatom_only,enl=enl) 382 else 383 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir_nonlop,lambda,& 384 & mpi_enreg,1,nnlout,paw_opt,signs,scwavef_out,tim_nonlop,cwavef,cwavef_out,& 385 & iatom_only=iatom_only) 386 end if 387 388 ! Post-processing if nonlop is called with specific options 389 if (signs==2) then 390 if (paw_opt<3) then 391 call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,cwavef_out,& 392 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 393 else 394 call dotprod_g(argr,argi,istwf_k,npw_k,cplex,cwavef,scwavef_out,& 395 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 396 end if 397 enlout(inlout)=argr 398 end if 399 if (signs==1.and.choice==1) then 400 call dotprod_g(argr,argi,istwf_k,npw_k,1,cwavef,cwavef,& 401 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 402 enlout(:)=enlout(:)+argr 403 end if 404 405 ! Write recognizable statements in log file 406 if (idtset<ndtset_test) then 407 write(std_out,'(a,i3,es24.16)') "TESTDFPT_df: ",iband,enlout(inlout) 408 else 409 write(std_out,'(a,i3,es24.16)') "TESTDFPT_dfpt:",iband,enlout(inlout) 410 end if 411 412 end if 413 414 end do ! End of loop on block of bands 415 416 ! Increment indexes (not used here because only one spin/kpoint) 417 icg=icg+npw_k*my_nspinor*nband_k 418 ikg=ikg+npw_k 419 420 end if ! Not my spin/kpoint 421 422 bdtot_index=bdtot_index+nband_k 423 424 !Memory deallocations 425 ABI_FREE(enl) 426 ABI_FREE(enlout) 427 ABI_FREE(lambda) 428 ABI_FREE(ph3d) 429 ABI_FREE(ffnl) 430 ABI_FREE(cwavef) 431 ABI_FREE(cwavef_out) 432 ABI_FREE(scwavef_out) 433 ABI_FREE(kg_k) 434 ABI_FREE(kpg_k) 435 ABI_FREE(ylm_k) 436 ABI_FREE(ylmgr_k) 437 ABI_FREE(ylm) 438 ABI_FREE(ylmgr) 439 call gs_hamk%free() 440 441 end subroutine nonlop_test