TABLE OF CONTENTS


ABINIT/m_nonlop_test [ Modules ]

[ Top ] [ Modules ]

NAME

  m_nonlop_test

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_nonlop_test
27 
28  implicit none
29 
30  private

ABINIT/nonlop_test [ Functions ]

[ Top ] [ 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=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

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